smfl.c

The Free Finite Element Package is a library which contains numerical methods required when working

      for (i = 0; i < mesh->number_of_triangles; i++)
	{
	  free(mesh->son[i]);
	}
      free(mesh->son);
    }

  free(mesh->boundary_edges);

  free (mesh);
}

/**
 * \~german 
 *  \brief Berechnet die Tranformationen vom Dreiecke der Triangulierung auf 
 *   das Refernezdreieck 
 *  
 *  \param points enthaelt die Informationen ueber die Punkte der Triangulierung
 *  in einem Format wie es von der Funktion smfl_read_point_data erzeugt wird.
 *
 *  \param triangles enthaelt die Informationen ueber die Knoten der Dreiecke 
 *  der Triangulierung in einem Format wie es von der 
 *  Funktion smfl_read_triangle_data erzeugt wird.
 *
*/
TRANSFORMATION *smfl_compute_backtransformations (const TRANSFORMATION *
						  trans,
						  const MEML_INT
						  number_of_triangles)
{
  TRANSFORMATION *transformation;
  MEML_INT i;
  MEML_FLOAT det;

  transformation =
    (TRANSFORMATION *) malloc (number_of_triangles * sizeof (TRANSFORMATION));

  for (i = 0; i < number_of_triangles; i++)
    {
      det = trans[i].det;

      transformation[i].B[0][0] = trans[i].B[1][1];
      transformation[i].B[1][0] = -trans[i].B[1][0];
      transformation[i].B[0][1] = -trans[i].B[0][1];
      transformation[i].B[1][1] = trans[i].B[0][0];

      transformation[i].det = 1 / det;
    }

  return (transformation);
}


/**
 * \~german 
 *  \brief Berechnet die Tranformationen vom Refernezdreieck auf die 
 *   Dreiecke der Triangulierung
 *  
 *  \param points enthaelt die Informationen ueber die Punkte der Triangulierung
 *  in einem Format wie es von der Funktion smfl_read_point_data erzeugt wird.
 *
 *  \param triangles enthaelt die Informationen ueber die Knoten der Dreiecke 
 *  der Triangulierung in einem Format wie es von der 
 *  Funktion smfl_read_triangle_data erzeugt wird.
 *
*/
TRANSFORMATION *smfl_compute_transformations (const VECTOR * points,
					      const INDEXARRAY * triangles)
{
  MEML_INT i;
  TRIANGLE the_triangle;
  T_POINTS the_points;
  TRANSFORMATION *transformation;
  MEML_INT *triangle_node = triangles->data;
  MEML_FLOAT *point = points->data;
  MEML_INT number_of_triangles;

  number_of_triangles = (triangles->dim) / 3;

  transformation =
    (TRANSFORMATION *) malloc (number_of_triangles * sizeof (TRANSFORMATION));

  for (i = 0; i < number_of_triangles; i++)
    {

      the_triangle = smfl_get_triangle (i, triangle_node);

      the_points = smfl_get_t_points (the_triangle, point);

      /*Bilde die Abbildung vom Referenzdreieck 
         auf das Dreieck */

      /* b00 b01
         b10 b11 */

      /* achtung habe hier am 24 die indiez vertauscht... */
      transformation[i].B[0][0] = the_points.a2[0] - the_points.a1[0];
      transformation[i].B[1][0] = the_points.a2[1] - the_points.a1[1];

      transformation[i].B[0][1] = the_points.a3[0] - the_points.a1[0];
      transformation[i].B[1][1] = the_points.a3[1] - the_points.a1[1];

      transformation[i].det =
	transformation[i].B[0][0] * transformation[i].B[1][1] -
	transformation[i].B[1][0] * transformation[i].B[0][1];

    }

  return (transformation);
}


/**
 * \~german 
 *  \brief lie遲 die Koordinaten der Gitterpunkte aus einer Datei ein
 *
 *   Die Funktion wurde geschrieben um eine einfachen Import von 
 *   Gitterdaten zu ermoeglichen.<br>
 *   Die Datei muss wie folgt aufgebaut sein<br><br>
 *   x1,y1<br>
 *   x2,y2<br>
 *    ... <br>
 *   xn,yn<br>
 *   <br><br>
 *  Leerzeilen oder aehnliche Abweichungen sind nicht zulae遡g.
 *  
 *  \param filename[] Name, bzw. Pfad+Name der Datei aus der die 
 *  Informationen gelesen werden sollen.
*/
VECTOR *smfl_read_point_data (char filename[])
{
  FILE *CSV;
  MEML_FLOAT testx, testy;
  MEML_INT counter = 0;
  MEML_INT number_of_points;
  VECTOR *rvector;

  number_of_points = 0;

  /*  CSV=fopen("POINT.data", "r"); */
  CSV = fopen (filename, "r");

  if (CSV == NULL)
    {
      fprintf (stderr, "Could not open %s for reading!\n", filename);
      exit (EXIT_FAILURE);
    }

  while ((fscanf (CSV, "%lf,%lf", &testx, &testy)) != EOF)
    number_of_points++;

  fclose (CSV);

  rvector = meml_vector_new (2 * number_of_points);

  CSV = fopen (filename, "r");

  while ((fscanf (CSV, "%lf,%lf", &testx, &testy)) != EOF)
    {
      rvector->data[(2 * counter)] = testx;
      rvector->data[(2 * counter) + 1] = testy;
      counter++;
    }

  rvector->dim = (number_of_points * 2);

  return (rvector);
}

/**
 * \~german 
 *  \brief lie遲 die Nummern der Gitterpunkte die auf dem Rand 
 *   des Gitters liegen aus einer Datei ein
 *
 *   Die Funktion wurde geschrieben um eine einfachen Import von 
 *   Gitterdaten zu ermoeglichen.<br>
 *   Die Datei muss wie folgt aufgebaut sein<br><br>
 *   b1<br>
 *   b2<br>
 *   ... <br>
 *   bn<br>
 *   <br><br>
 *  Leerzeilen oder aehnliche Abweichungen sind nicht zulae遡g.
 *  
 *  \param filename[] Name, bzw. Pfad+Name der Datei aus der die 
 *  Informationen gelesen werden sollen.
*/
INDEXARRAY *smfl_read_boundary_data (char filename[])
{
  /*  (MEML_INT ** b_points, */
  FILE *CSV;
  long int node1;
  MEML_INT counter = 0;
  INDEXARRAY *rvector;
  MEML_INT number_of_b_points = 0;

  /* Die Dreicke werden in einem array of int abgelegt 
     Speicher fuer den Array wird hier besorgt 
     T=[ 2 3 13 17 8 9 .. ]
     2 3 13 sind die Punkte des ersten Dreiecks
     17 8 9 sind die Punkte des folgenden Dreiecks
   */

  CSV = fopen (filename, "r");

  if (CSV == NULL)
    {
      fprintf (stderr, "Could not open %s for reading!\n", filename);
      exit (EXIT_FAILURE);
    }

  while ((fscanf (CSV, "%ld", &node1)) != EOF)
    number_of_b_points++;

  fclose (CSV);

  rvector = meml_indexarray_new (number_of_b_points);

  CSV = fopen (filename, "r");

  while ((fscanf (CSV, "%ld", &node1)) != EOF)
    {
      rvector->data[counter] = node1;
      counter++;
    }

  rvector->dim = number_of_b_points;

  return (rvector);

}


MESH * smfl_read_gts_meshdata(char kind, char filename[])
{
  FILE * CSV;
  MEML_FLOAT testx, testy, foo;
  char * trash1, *trash2, *trash3, *trash4;
  VECTOR *points;
  int NumberOfPoints=0, NumberOfTriangles=0, NumberOfEgdes=0;
  int i,j,k, ** edges;
  int triangletemp[3];
  INDEXARRAY *triangles;
  MESH * mesh;

  CSV = fopen (filename, "r");

  if (CSV == NULL)
    {
      fprintf (stderr, "Could not open %s for reading!\n", filename);
      exit (EXIT_FAILURE);
    }

  trash1 = (char*) calloc(100,sizeof(char));
  trash2 = (char*) calloc(100,sizeof(char));
  trash3 = (char*) calloc(100,sizeof(char));
  trash4 = (char*) calloc(100,sizeof(char));

  fscanf (CSV, "%d %d %d  %s %s %s %s", 
	  &NumberOfPoints, &NumberOfEgdes, &NumberOfTriangles, 
	  trash1, trash2, trash3, trash4);

  free(trash1);  free(trash2);  free(trash3);  free(trash4);

  printf("NumberOfPoints %d NumberOfEgdes %d NumberOfTriangles %d \n", 
	 NumberOfPoints, NumberOfEgdes, NumberOfTriangles);

  points = meml_vector_new(NumberOfPoints*2);

  for (i=0;i<NumberOfPoints;i++)
    {
      fscanf (CSV, "%lf %lf %lf", &testx, &testy, &foo);
      points->data[2*i] = testx; 
      points->data[2*i +1] = testy;
    }

  meml_vector_print(points);

  edges = (int **) calloc(NumberOfEgdes,sizeof(int*));

  for (i=0;i<NumberOfEgdes;i++)
    {
      edges[i] = (int *) calloc(2,sizeof(int) );
      fscanf (CSV, "%d %d", &(edges[i][0]), &(edges[i][1]));
    }

  triangles = meml_indexarray_new(NumberOfTriangles*3);

  for (i=0;i<NumberOfTriangles;i++)
    {
      fscanf (CSV, "%d %d %d", &(triangletemp[0]),&(triangletemp[1]),&(triangletemp[2]));

      triangles->data[3*i] = edges[triangletemp[0]-1][0];
      triangles->data[3*i+1]  = edges[triangletemp[0]-1][1];
      triangles->data[3*i+2]  = -1;
      
      for (j=1;j<3;j++)
      {
	  for (k=0;k<2;k++)
	      if ( ( edges[triangletemp[j]-1][k] != triangles->data[3*i] ) &&  
		   ( edges[triangletemp[j]-1][k] != triangles->data[3*i] ) )
	      {
		  triangles->data[3*i+2]  = edges[triangletemp[j]-1][k];
	      }
      }
      
      if (triangles->data[3*i+2]  == -1 )
      {
	  fprintf(stderr,"unkown error in smfl_read_gts_meshdata!\n"
		  "Maybe the gts file is corrupt!\n");
	  exit(EXIT_FAILURE);
      }
    }
  
  for (i=0;i<NumberOfEgdes;i++)
    {
      free(edges[i]);
    }
  free(edges);

  fclose (CSV);

  mesh = smfl_mesh_new (kind, NULL, points,triangles);

  meml_vector_free(points);
  meml_indexarray_free(triangles);

  return(mesh);
}


MESH * smfl_read_msh_meshdata(char kind, char filename[])
{
  FILE * CSV;
  MEML_FLOAT testx, testy, foo;
  int foo1, foo2, foo3, elementNo, elementtype;
  char * trash1, *trash2, *trash3, *trash4;
  VECTOR *points;
  int NumberOfPoints=0, NumberOfTriangles=0;
  int i,j, currTriangle=0;
  int triangletemp[3];
  INDEXARRAY *triangles, *renumber,*trianglesTemp;
  MESH * mesh;

  CSV = fopen (filename, "r");

  if (CSV == NULL)
    {
      fprintf (stderr, "Could not open %s for reading!\n", filename);
      exit (EXIT_FAILURE);
    }

  trash1 = (char*) calloc(100,sizeof(char));
  trash2 = (char*) calloc(100,sizeof(char));
  trash3 = (char*) calloc(100,sizeof(char));
  trash4 = (char*) calloc(100,sizeof(char));

  fscanf (CSV, "%s",trash1); 
  fscanf (CSV, "%d",&NumberOfPoints); 

  points = meml_vector_new(NumberOfPoints*2);
  renumber = meml_indexarray_new(NumberOfPoints);

  for (i=0;i<NumberOfPoints;i++)
    {
      fscanf (CSV, "%d %lf %lf %lf", &elementNo, &testx, &testy, &foo);
      points->data[2*i] = testx; 
      points->data[2*i +1] = testy;
      renumber->data[i]=elementNo;
    }

  fscanf (CSV, "%s",trash1); 
  fscanf (CSV, "%s",trash1); 
  fscanf (CSV, "%d",&NumberOfTriangles); 

  trianglesTemp = meml_indexarray_new(NumberOfTriangles*3);

  for (i=0;i<trianglesTemp->dim/3;i++)
    {
	fscanf (CSV, "%d %d %d %d %d ", &elementNo, 
		&elementtype, &foo1, &foo2, &foo3);

	if (elementtype == 1) 
	  fscanf (CSV, "%d %d \n", &foo1, &foo2);

	if (elementtype == 2)
	{
	  fscanf (CSV, "%d %d %d \n", 
		  &(triangletemp[0]),&(triangletemp[1]),&(triangletemp[2]));
	  
	  for (j=0;j<renumber->dim;j++)
	    {
	      if (renumber->data[j] == triangletemp[0])
		{
		  trianglesTemp->data[3*currTriangle] = j+1;		    
		}
	      if (renumber->data[j] == triangletemp[1])
		{
		  trianglesTemp->data[3*currTriangle+1] = j+1;		    
		}
	      if (renumber->data[j] == triangletemp[2])
		{
		  trianglesTemp->data[3*currTriangle+2] = j+1;		    
		}
	    }
	  currTriangle++;
	}
	else
	{
	    NumberOfTriangles--;
	}
	

    }
  fclose (CSV);

  meml_indexarray_free(renumber);

  triangles = meml_indexarray_new(NumberOfTriangles*3);

  for (i=0;i<triangles->dim;i++)
  {
      triangles->data[i] = trianglesTemp->data[i];
  }
  
  meml_indexarray_free(trianglesTemp);

  mesh = smfl_mesh_new (kind, NULL, points,triangles);

  meml_vector_free(points);
  meml_indexarray_free(triangles);

  return(mesh);
}


/**
 * \~german 
 *  \brief lie遲 Informationen ueber den Aufbau der einzelnen 
 *   Dreiecke aus einer Datei ein.
 *
 *   Die Funktion wurde geschrieben um eine einfachen Import von 
 *   Gitterdaten zu ermoeglichen.<br><br>
 *   Ein Dreieck wird durch die Nummer seiner drei Eckpunkte 
 *   definiert.<br>
 *   Darueberhinaus exportieren viele Programme noch Informationen 
 *   ueber unterschiedliche Layer, die jedoch von dieser Funktion 
 *   nicht ausgewertet werden.<br>
 *   <br>
 *   Die Datei muss wie folgt aufgebaut sein<br><br>
 *   \f$T1_1,T1_2,T1_3,layer\f$<br>
 *   \f$T2_1,T2_2,T2_3,layer\f$<br>
 *   ... <br>
 *   \f$Tn_1,Tn_2,Tn_3,layer\f$<br>
 *   <br><br>
 *  Leerzeilen oder aehnliche Abweichungen sind nicht zulae遡g.<br>
 *  Die Information layer wird verworfen
 *  
 *  \param filename[] Name, bzw. Pfad+Name der Datei aus der die 
 *  Informationen gelesen werden sollen.
*/
INDEXARRAY *smfl_read_triangle_data (char filename[])
{
  FILE *CSV;
  long int node1, node2, node3, foo;
  MEML_INT counter = 0;
  MEML_INT number_of_points = 0;
  INDEXARRAY *rvector;

  /* Die Dreicke werden in einem array of int abgelegt 
     Speicher fuer den Array wird hier besorgt 
     T=[ 2 3 13 17 8 9 .. ]
     2 3 13 sind die Punkte des ersten Dreiecks
     17 8 9 sind die Punkte des folgenden Dreiecks
   */

  CSV = fopen (filename, "r");

  if (CSV == NULL)
    {
      fprintf (stderr, "Could not open %s for reading!\n", filename);
      exit (EXIT_FAILURE);
    }

  while ((fscanf (CSV, "%ld,%ld, %ld, %ld",
		  &node1, &node2, &node3, &foo)) != EOF)
    number_of_points = number_of_points + 3;

  fclose (CSV);

  rvector = meml_indexarray_new (number_of_points);

  CSV = fopen (filename, "r");

  while ((fscanf (CSV, "%ld,%ld, %ld, %ld",
		  &node1, &node2, &node3, &foo)) != EOF)
    {
      rvector->data[counter] = node1;
      rvector->data[counter + 1] = node2;
      rvector->data[counter + 2] = node3;
      counter = counter + 3;
    }

  fclose (CSV);

  /* Achtung das Bedeutet noch aerger 
     (*number_of_points)=(*number_of_points)/3; */

  return (rvector);

}