smfl.c

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

	    }
	}
      /* If only the longest edge is divided, 
         its midpoint is joined with the opposing corner. */
      else if (number_of_mark_edges == 1)
	{
	  /* gruen */
	  triangles->data[zeiger * 3] = kanten[welche][0];
	  triangles->data[zeiger * 3 + 1] = gegenueber;
	  triangles->data[zeiger * 3 + 2] = mittelpunkte[welche];
	  mesh->son[i][0] = zeiger;
	  zeiger += 1;

	  triangles->data[zeiger * 3] = kanten[welche][1];
	  triangles->data[zeiger * 3 + 1] = gegenueber;
	  triangles->data[zeiger * 3 + 2] = mittelpunkte[welche];
	  mesh->son[i][1] = zeiger;
	  mesh->son[i][2] = -1;
	  zeiger += 1;
	}
    }



  newMesh = smfl_mesh_new (mesh->kind, boundary, points, triangles);

  if (mesh->kind == 'l')
  {
      newMesh->restriction = meml_matrix_convert (CS, restriction);
  }
  else
  {
      meml_matrix_free (restriction);

      restriction =
	  meml_matrix_new (LS,  
			   mesh->number_of_points+mesh->number_of_adof,
			   newMesh->number_of_points+newMesh->number_of_adof);

      /* alte Punkte uebernehmen */
      for (i = 0; i < mesh->number_of_points; i++)
      {
	  meml_matrix_element_set_f (restriction,i,i, 1);
      }

      /* neuen Teil durchgehen */
      zeiger = mesh->points->dim / 2 + 1;

      for (i = 0; i < mark_edges->dim; i++)
      {
	  /* alter adof ist neuer Punkt, also permutieren */
	  if (mark_edges->data[i] > 0)
	  {
	      meml_matrix_element_set_f (restriction, mesh->number_of_points + i,zeiger - 1 ,1);
	      zeiger++;
	  }
      }

      for (i = 0; i < newMesh->number_of_edges; i++)
      {
	  /* gab es die kante schon frueher ? */
	  if ( (newMesh->edges[i][0] <= mesh->number_of_points) &&
	       (newMesh->edges[i][1] <= mesh->number_of_points) )
	  {

	      j=smfl_in_liste_fast (newMesh->edges[i][0],newMesh->edges[i][1],
				    mesh->edges, 
				    mesh->number_of_edges);

	      meml_matrix_element_set_f (restriction, 
					 mesh->number_of_points + j,
					 newMesh->number_of_points + i,
					 1);
	  }
      }



      for (i = 0; i < newMesh->number_of_edges; i++)
      {
	  if ( (newMesh->edges[i][0] > mesh->number_of_points) ||
	       (newMesh->edges[i][1] > mesh->number_of_points) )
	  {
	      MEML_FLOAT * values;
	      MEML_INT * rows;
	      MEML_INT entries,where,where2;

	      /* wenn ein neuer Punkt drin ist : nein */
	      /* dann aus den beiden punkten interpolieren */
	      /* dazu alte adof-nummer ermitteln */
	      if (newMesh->edges[i][0] <= mesh->number_of_points)
	      {
		  meml_matrix_element_set_f (restriction, 
					     newMesh->edges[i][0]-1 ,
					     newMesh->number_of_points + i, 
					     0.5);
		  ssml_ls_col_get (restriction->data.flist_sparse_matrix,newMesh->edges[i][1]-1,
				   &values, &rows,&entries);

		  for (k=0;k<entries;k++)
		  {
		      if (values[k]==1)
		      {
			  where =rows[k]; 		  
		      }
		  }
		  meml_matrix_element_set_f (restriction,where,newMesh->number_of_points + i, 0.5);
		  free(values);
		  free(rows);
	      }
	      else if (newMesh->edges[i][1] <= mesh->number_of_points) 
	      {
		  meml_matrix_element_set_f (restriction, 
					     newMesh->edges[i][1]-1 ,
					     newMesh->number_of_points + i, 
					     0.5);

		  ssml_ls_col_get (restriction->data.flist_sparse_matrix,newMesh->edges[i][0]-1,
				   &values, &rows,&entries);

		  for (k=0;k<entries;k++)
		  {
		      if (values[k]==1)
		      {
			  where =rows[k]; 		  
		      }
		  }
		  meml_matrix_element_set_f (restriction, where, newMesh->number_of_points + i, 0.5);
		  free(values);
		  free(rows);
	      }
	      else
	      {
		  ssml_ls_col_get (restriction->data.flist_sparse_matrix,newMesh->edges[i][0]-1,
				   &values, &rows,&entries);
		  
		  for (k=0;k<entries;k++)
		  {
		      if (values[k]==1)
		      {
			  where =rows[k]; 		  
		      }
		  }
		  free(values);
		  free(rows);		  
		  ssml_ls_col_get (restriction->data.flist_sparse_matrix,newMesh->edges[i][1]-1,
				   &values, &rows,&entries);
		  
		  for (k=0;k<entries;k++)
		  {
		      if (values[k]==1)
		      {
			  where2 =rows[k]; 		  
		      }
		  }
		  free(values);
		  free(rows);		  

		  meml_matrix_element_set_f (restriction, where,newMesh->number_of_points + i, 0.5);
		  meml_matrix_element_set_f (restriction, where2,newMesh->number_of_points + i, 0.5);
	      }
	  }
      }
      newMesh->restriction = meml_matrix_convert (CS, restriction);
  }

  meml_matrix_free (restriction);
  meml_vector_free(mark_edges);
  meml_vector_free (points);
  meml_indexarray_free (boundary);
  meml_indexarray_free (triangles);

  return (newMesh);
}

MESH *  smfl_mesh_global_refine(MESH * mesh)
{
  INDEXARRAY * torefine;
  MEML_INT i;
  MESH * newMesh;

  torefine=meml_indexarray_new(mesh->number_of_triangles);

  for (i=0;i<mesh->number_of_triangles;i++)
    torefine->data[i]=i;

  newMesh=smfl_mesh_refine (mesh,torefine);

  meml_indexarray_free(torefine);

  return(newMesh);
}


void smfl_mesh_guess_normals (MESH * mesh)
{
  MEML_INT i, p1, p2;
  MEML_FLOAT x1, y1, x2, y2, x3, y3, n1, n2, norm;
  MEML_INT which_edge, which_triangle;
  TRIANGLE dreieck;

  mesh->normal_vector =
    (VECTOR **) malloc (mesh->number_of_boundary_edges * sizeof (VECTOR *));

  for (i = 0; i < mesh->number_of_boundary_edges; i++)
    mesh->normal_vector[i] = meml_vector_new (2);

  for (i = 0; i < mesh->number_of_boundary_edges; i++)
    {
      /* remember it the triangle = [1 5 3] is mapped on the 
         reference triangle 
         1 is (0,0) ; 5 is (1,0) and 3 is (0,1) */

      which_edge = mesh->boundary_edges[i];

      p1 = mesh->edges[which_edge][0];
      p2 = mesh->edges[which_edge][1];

      x1 = mesh->points->data[2 * (p1 - 1)];
      y1 = mesh->points->data[2 * (p1 - 1) + 1];

      x2 = mesh->points->data[2 * (p2 - 1)];
      y2 = mesh->points->data[2 * (p2 - 1) + 1];

      n1 = (y2 - y1);
      n2 = -(x2 - x1);

      norm = sqrt ((y2 - y1) * (y2 - y1) + (x2 - x1) * (x2 - x1));

      n1 = n1 / norm;
      n2 = n2 / norm;

      /* a boundary edge has only one triangle */
      which_triangle = mesh->edge2triangle[which_edge][0];
      dreieck = smfl_get_triangle (which_triangle, mesh->triangles->data);

      /* add the normal vector to one of the edge points */

      x2 = x1 - n1;
      y2 = y1 - n2;

      x1 = x1 + n1;
      y1 = y1 + n2;

      /* get the 3 point of the triangle */
      if ((dreieck.p[0] != p1) && (dreieck.p[0] != p2))
	{
	  x3 = mesh->points->data[2 * (dreieck.p[0] - 1)];
	  y3 = mesh->points->data[2 * (dreieck.p[0] - 1) + 1];
	}
      else if ((dreieck.p[1] != p1) && (dreieck.p[1] != p2))
	{
	  x3 = mesh->points->data[2 * (dreieck.p[1] - 1)];
	  y3 = mesh->points->data[2 * (dreieck.p[1] - 1) + 1];
	}
      else if ((dreieck.p[2] != p1) && (dreieck.p[2] != p2))
	{
	  x3 = mesh->points->data[2 * (dreieck.p[2] - 1)];
	  y3 = mesh->points->data[2 * (dreieck.p[2] - 1) + 1];
	}
      else
	{
	  x3 = 0;
	  y3 = 0;
	  fprintf (stderr, "Error in smfl_mesh_new : "
		   "One or more triangles are ill-defined!\n");
	  exit (EXIT_FAILURE);
	}

      /* does n have the correct direction ? */
      /* if the correct n will have the greater distance to 
         the 3 point of the triangle */
      if (sqrt ((x3 - x1) * (x3 - x1) + (y3 - y1) * (y3 - y1)) <
	  sqrt ((x3 - x2) * (x3 - x2) + (y3 - y2) * (y3 - y2)))
	{
	  n1 = -n1;
	  n2 = -n2;
	}

      mesh->normal_vector[i]->data[0] = n1;
      mesh->normal_vector[i]->data[1] = n2;

    }
}

static void h_calculating(MESH * mesh)
{
  MEML_INT i, number_of_triangles = mesh->number_of_triangles;
  VECTOR * h;
  MEML_FLOAT temp[2];
  TRANSFORMATION *transformation = mesh->transformation;

  mesh->ht = (MEML_FLOAT *) calloc (number_of_triangles, sizeof(MEML_FLOAT));
  mesh->hmax = 0;
  /* um inf zu erzeugen */
  mesh->hmin = 1.0/0.0;

  h = meml_vector_new(3);

  for (i = 0;i<number_of_triangles;i++)
    {
      temp[0] = (transformation[i].B[0][0] - transformation[i].B[0][1]);
      temp[0] = temp[0]*temp[0];

      temp[1] = (transformation[i].B[1][0] - transformation[i].B[1][1]);
      temp[1] = temp[1]*temp[1];

      h->data[2] = sqrt(temp[0]+temp[1]);

      h->data[0] = sqrt(transformation[i].B[0][0]*transformation[i].B[0][0] + 
		  transformation[i].B[1][0] *transformation[i].B[1][0] );

      h->data[1] = sqrt(transformation[i].B[0][1]*transformation[i].B[0][1] + 
		  transformation[i].B[1][1] *transformation[i].B[1][1] );
      
      mesh->ht[i] = meml_vector_norm_max(h);
      
      if ( mesh->ht[i] > mesh->hmax)
	{
	  mesh->hmax = mesh->ht[i];
	}
      if ( mesh->ht[i] < mesh->hmin)
	{
	  mesh->hmin = mesh->ht[i];
	}

    }

  meml_vector_free(h);
}

void smfl_boundary_find(MESH *mesh)
{
    INDEXARRAY * temp;
    MEML_INT i,welcheKante,wieviele=0,zaehler=0;
    
    temp =  meml_indexarray_new (mesh->number_of_points);
    
    for (i=0;i<mesh->number_of_boundary_edges;i++)
    {
	welcheKante=mesh->boundary_edges[i];
	
	temp->data[mesh->edges[welcheKante][0]-1] = 1;
	temp->data[mesh->edges[welcheKante][1]-1] = 1;
    }
	
    for (i=0;i<temp->dim;i++)
    {
	if (temp->data[i] == 1)
	    wieviele++;
    }

    mesh->boundary = meml_indexarray_new(wieviele);

    for (i = 0; i < temp->dim; i++)
      {
	  if (temp->data[i] == 1)
	  {
	      mesh->boundary->data[zaehler]=i;
	      zaehler++;
	  }
      }

    meml_indexarray_free(temp);
}

/**
 * \brief Creates a new mesh based on the user data.
 *
 * \param boundary the number of the nodes which belongs to the boundary
 * \param points the xy coordinates of the mesh nodes 
 * \param triangles information about the triangulation 
 *
 */
MESH *smfl_mesh_new (char kind, INDEXARRAY * boundary, VECTOR * points,
		     INDEXARRAY * triangles)
{
  MESH *new_mesh;

  new_mesh = (MESH *) calloc (1, sizeof (MESH));

  new_mesh->kind = kind;
  new_mesh->restriction = NULL;

  if (boundary != NULL)
  {
      new_mesh->boundary = meml_indexarray_copy (boundary);
      new_mesh->number_of_boundary_points = boundary->dim;
  }
  else
  {
      new_mesh->boundary = NULL;
  }

  new_mesh->points = meml_vector_copy (points);
  new_mesh->triangles = meml_indexarray_copy (triangles);

  new_mesh->number_of_points = points->dim / 2;
  new_mesh->number_of_triangles = triangles->dim / 3;
  new_mesh->number_of_edges = new_mesh->number_of_triangles * 3;
  /*
    new_mesh->number_of_triangles + new_mesh->number_of_points - 1;
  */

  /* compute transformations */
  new_mesh->transformation = smfl_compute_transformations (points, triangles);

  new_mesh->backtransformation =
    smfl_compute_backtransformations (new_mesh->transformation,
				      new_mesh->number_of_triangles);

  triangles2edges (new_mesh);

  new_mesh->normal_vector = NULL;
  new_mesh->number_of_adof = 0;

  /* bubbles, local quadratic or linear ? */
  switch (kind)
    {
    case 'b':
      new_mesh->degrees_of_freedom_per_triangle = 4;
      break;
    case 'q':
      new_mesh->degrees_of_freedom_per_triangle = 6;
      /* calculating triangles_adof and additional_degrees_of_freedom */
      calculating_triangles_adof (new_mesh);
      break;
    default:
      new_mesh->degrees_of_freedom_per_triangle = 3;
    }

  /* calculating point2triangle */
  calculating_points2triangle (new_mesh);

  if (new_mesh->boundary == NULL)
  {
      smfl_boundary_find(new_mesh);
      new_mesh->number_of_boundary_points = new_mesh->boundary->dim;
  }

  /* Information ueber die Nachbarpunkte eines Knoten */
  getting_neighbours (new_mesh);

  h_calculating(new_mesh);

  new_mesh->son=NULL;

  return (new_mesh);
}

/**
 * \brief Frees the memory used by mesh
 *
 */
void smfl_mesh_free (MESH * mesh)
{
  MEML_INT i;

  for (i = 0; i < mesh->number_of_points; i++)
    {
      meml_indexarray_free (mesh->point2triangle[i]);
    }
  free (mesh->point2triangle);

  for (i = 0; i < mesh->number_of_edges; i++)
    {
      free (mesh->edge2triangle[i]);
      free (mesh->edges[i]);
    }

  free (mesh->edges);
  free (mesh->edge2triangle);

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

  if (mesh->triangles_adof != NULL)
    for (i = 0; i < mesh->number_of_triangles; i++)
      {
	free (mesh->triangles_adof[i]);
      }

  free (mesh->transformation);
  free (mesh->backtransformation);
  free (mesh->triangle2edge);
  free (mesh->triangles_adof);
  free (mesh->ht);

  meml_vector_free (mesh->points);
  if (mesh->additional_degrees_of_freedom != NULL)
    meml_vector_free (mesh->additional_degrees_of_freedom);

  meml_indexarray_free (mesh->boundary);
  meml_indexarray_free (mesh->triangles);

  for (i = 0; i < mesh->number_of_points; i++)
    meml_indexarray_free(mesh->neighbours[i]);

  free (mesh->neighbours);
  free (mesh->number_of_neighbours);

  if (mesh->restriction != NULL)
    meml_matrix_free (mesh->restriction);

  if (mesh->normal_vector !=NULL)
    {
      for (i = 0; i < mesh->number_of_boundary_edges; i++)
	meml_vector_free(mesh->normal_vector[i]);
      free(mesh->normal_vector);
    }

  if (mesh->son !=NULL)
    {