smfl.c

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

{
  MEML_INT i;
  MEML_INT position;

  new_mesh->number_of_neighbours =
    (MEML_INT *) calloc (new_mesh->number_of_points, sizeof (MEML_INT));

  new_mesh->neighbours =
    (INDEXARRAY **) calloc (new_mesh->number_of_points,
			    sizeof (INDEXARRAY *));

  for (i = 0; i < new_mesh->number_of_edges; i++)
    {
      new_mesh->number_of_neighbours[new_mesh->edges[i][0] - 1]++;
      new_mesh->number_of_neighbours[new_mesh->edges[i][1] - 1]++;
    }

  for (i = 0; i < new_mesh->number_of_points; i++)
    new_mesh->neighbours[i] =
      meml_indexarray_new (new_mesh->number_of_neighbours[i]);

  for (i = 0; i < new_mesh->number_of_edges; i++)
    {
      position = 0;
      while (new_mesh->neighbours[new_mesh->edges[i][0] - 1]->
	     data[position] != 0)
	position++;
      new_mesh->neighbours[new_mesh->edges[i][0] - 1]->data[position] =
	new_mesh->edges[i][1];

      position = 0;
      while (new_mesh->neighbours[new_mesh->edges[i][1] - 1]->
	     data[position] != 0)
	position++;
      new_mesh->neighbours[new_mesh->edges[i][1] - 1]->data[position] =
	new_mesh->edges[i][0];
    }
}

static MESH *smfl_mesh_copy (const MESH * mesh)
{
  MESH *meshCopy;
  MEML_INT i, j;

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

  meshCopy->kind = mesh->kind;

  meshCopy->number_of_points = mesh->number_of_points;
  meshCopy->number_of_triangles = mesh->number_of_triangles;
  meshCopy->number_of_boundary_points = mesh->number_of_boundary_points;
  meshCopy->number_of_edges = mesh->number_of_edges;
  meshCopy->number_of_boundary_edges = mesh->number_of_boundary_edges;
  meshCopy->number_of_adof = mesh->number_of_adof;
  meshCopy->degrees_of_freedom_per_triangle =
    mesh->degrees_of_freedom_per_triangle;

  meshCopy->boundary_edges =
    (MEML_INT *) calloc (meshCopy->number_of_boundary_edges,
			 sizeof (MEML_INT));

  for (i = 0; i < meshCopy->number_of_boundary_edges; i++)
    meshCopy->boundary_edges[i] = mesh->boundary_edges[i];

  meshCopy->number_of_neighbours =
    (MEML_INT *) calloc (meshCopy->number_of_points, sizeof (MEML_INT));
  meshCopy->normal_vector =
    (VECTOR **) calloc (meshCopy->number_of_points, sizeof (VECTOR *));

  for (i = 0; i < meshCopy->number_of_points; i++)
    {
      meshCopy->number_of_neighbours[i] = mesh->number_of_neighbours[i];
      meshCopy->normal_vector[i] = meml_vector_copy (mesh->normal_vector[i]);
    }

  meshCopy->points = meml_vector_copy (mesh->points);
  meshCopy->additional_degrees_of_freedom =
    meml_vector_copy (mesh->additional_degrees_of_freedom);

  meshCopy->boundary = meml_indexarray_copy (mesh->boundary);
  meshCopy->triangles = meml_indexarray_copy (mesh->triangles);

  if (meshCopy->kind == 'q')
    {
      meshCopy->triangles_adof =
	(MEML_INT **) calloc (meshCopy->number_of_points,
			      sizeof (MEML_INT *));

      for (i = 0; i < meshCopy->number_of_triangles; i++)
	{
	  meshCopy->triangles_adof[i] =
	    (MEML_INT *) calloc (3, sizeof (MEML_INT));
	  for (j = 0; j < 3; j++)
	    {
	      meshCopy->triangles_adof[i][j] = mesh->triangles_adof[i][j];
	    }
	}
    }
  else
    meshCopy->triangles_adof = NULL;

  meshCopy->edges =
    (MEML_INT **) calloc (meshCopy->number_of_edges, sizeof (MEML_INT *));
  meshCopy->edge2triangle =
    (MEML_INT **) calloc (meshCopy->number_of_edges, sizeof (MEML_INT *));

  for (i = 0; i < meshCopy->number_of_edges; i++)
    {
      meshCopy->edges[i] = (MEML_INT *) calloc (2, sizeof (MEML_INT));
      meshCopy->edge2triangle[i] = (MEML_INT *) calloc (2, sizeof (MEML_INT));
      for (j = 0; j < 2; j++)
	{
	  meshCopy->edges[i][j] = mesh->edges[i][j];
	  meshCopy->edge2triangle[i][j] = mesh->edge2triangle[i][j];
	}
    }

  meshCopy->triangle2edge =
    (MEML_INT **) calloc (meshCopy->number_of_triangles, sizeof (MEML_INT *));

  for (i = 0; i < meshCopy->number_of_triangles; i++)
    {
      meshCopy->triangle2edge[i] = (MEML_INT *) calloc (3, sizeof (MEML_INT));
      for (j = 0; j < 3; j++)
	{
	  meshCopy->triangle2edge[i][j] = meshCopy->triangle2edge[i][j];
	}
    }

  meshCopy->point2triangle =
    (INDEXARRAY **) calloc (mesh->number_of_points, sizeof (INDEXARRAY *));
  meshCopy->neighbours =
    (INDEXARRAY **) calloc (mesh->number_of_points, sizeof (INDEXARRAY *));

  for (i = 0; i < mesh->number_of_points; i++)
    {
      meshCopy->point2triangle[i] =
	meml_indexarray_copy (mesh->point2triangle[i]);
      meshCopy->neighbours[i] = meml_indexarray_copy (mesh->neighbours[i]);
    }

  meshCopy->transformation =
    smfl_compute_transformations (meshCopy->points, meshCopy->triangles);

  if (meshCopy->restriction != NULL)
    meshCopy->restriction = meml_matrix_copy (mesh->restriction);

  return (meshCopy);
}


MESH *smfl_mesh_refine (MESH * mesh, const INDEXARRAY * torefine)
{
  VECTOR *mark_edges;
  MEML_INT i, j, k, kante, number_of_mark_edges = 0;
  ME_MATRIX *refine;
  MEML_INT newpoints = 0, newbpoints = 0, checkup = 0;
  register MEML_INT number_of_triangles = mesh->number_of_triangles;
  MEML_INT welche, wieviele, zeiger, bzeiger, zaehler;
  MEML_FLOAT laengste, laenge, x[2], y[2];
  VECTOR *points;
  INDEXARRAY *triangles, *boundary;
  ME_MATRIX *restriction;
  MEML_INT kanten[3][2], mittelpunkte[3], gegenueber = -1, t[3];
  MESH *newMesh;

  mesh->son =
    (MEML_INT **) calloc (mesh->number_of_triangles, sizeof (MEML_INT *));

  for (i = 0; i < mesh->number_of_triangles; i++)
    {
      mesh->son[i] = (MEML_INT *) calloc (4, sizeof (MEML_INT));
    }

  /* refine = meml_matrix_new (CS, mesh->number_of_triangles, 1); */
  mark_edges = meml_vector_new_zeros (mesh->number_of_edges);

  /* Either divide all edges of the selected triangles in half 
     (regular refinement) */

  for (i = 0; i < torefine->dim; i++)
    {
      /* meml_matrix_element_set_f (refine, torefine->data[i], 1, 1); */
      for (j = 0; j < 3; j++)
	{
	  kante = mesh->triangle2edge[torefine->data[i]][j];
	  mark_edges->data[kante] = 1;
	}
    }

  /* Divide the longest edge of any triangle that has a divided edge */
  /* Repeat step this until no further edges are divided */
  while (checkup == 0)
    {
      checkup = 1;

      for (i = 0; i < number_of_triangles; i++)
	{
	  /* wieviele kanten sind makiert */
	  number_of_mark_edges = 0;
	  for (j = 0; j < 3; j++)
	    {
	      kante = mesh->triangle2edge[i][j];
	      if (mark_edges->data[kante] == 1)
		number_of_mark_edges++;
	    }

	  /* wenn es nur 1 oder 2 sind, muss geprueft werden, 
	     ob die laengste dabei ist */
	  if ((number_of_mark_edges == 2) || (number_of_mark_edges == 1))
	    {
	      /* welche ist die laengste ? */
	      laengste = 0;
	      welche = -1;

	      for (j = 0; j < 3; j++)
		{
		  kante = mesh->triangle2edge[i][j];
		  x[0] = mesh->points->data[2 * (mesh->edges[kante][0] - 1)];
		  x[1] = mesh->points->data[2 * (mesh->edges[kante][1] - 1)];
		  y[0] =
		    mesh->points->data[2 * (mesh->edges[kante][0] - 1) + 1];
		  y[1] =
		    mesh->points->data[2 * (mesh->edges[kante][1] - 1) + 1];
		  laenge =
		    (x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] -
								     y[1]);

		  if (laenge > laengste)
		    {
		      laengste = laenge;
		      welche = j;
		    }
		}

	      kante = mesh->triangle2edge[i][welche];

	      /* wenn die laengste noch nicht makiert war, 
	         wird sie makiert und 
	         checkup wird zurueck gesetzt  */
	      if (mark_edges->data[kante] != 1)
		{
		  mark_edges->data[kante] = 1;
		  checkup = 0;
		}
	    }
	}
    }

  /* Introduce new points of all divided edges, 
     and replace all divided entries in
     mesh->boundary by two new entries. */

  for (i = 0; i < mark_edges->dim; i++)
    {
      if (mark_edges->data[i] == 1)
	{
	  newpoints++;
	  for (j = 0; j < mesh->number_of_boundary_edges; j++)
	    {
	      if (mesh->boundary_edges[j] == i)
		{
		  newbpoints++;
		  break;
		}
	    }
	}
    }

  /* loeschen  */

  /* copy old data */
  points = meml_vector_new (2 * (mesh->number_of_points + newpoints));
  boundary =
    meml_indexarray_new (mesh->boundary->dim + newbpoints);

  for (i = 0; i < mesh->points->dim; i++)
    points->data[i] = mesh->points->data[i];

  for (i = 0; i < mesh->boundary->dim; i++)
    boundary->data[i] = mesh->boundary->data[i];

  
  restriction =
      meml_matrix_new_eye (LS, mesh->number_of_points, points->dim / 2);
  
  /* add new data */

  /* loeschen  */

  zeiger = mesh->points->dim / 2 + 1;
  bzeiger = mesh->boundary->dim;
  for (i = 0; i < mark_edges->dim; i++)
    {
      if (mark_edges->data[i] == 1)
	{
	  /* zuordnung zum neuen Punkt fuer dreiecke speichern */
	  mark_edges->data[i] = zeiger;

	  /* mittelpunkt bestimmen */
	  x[0] = mesh->points->data[2 * (mesh->edges[i][0] - 1)];
	  x[1] = mesh->points->data[2 * (mesh->edges[i][1] - 1)];
	  y[0] = mesh->points->data[2 * (mesh->edges[i][0] - 1) + 1];
	  y[1] = mesh->points->data[2 * (mesh->edges[i][1] - 1) + 1];

	  points->data[2 * (zeiger - 1)] = (x[0] + x[1]) / 2;
	  points->data[2 * (zeiger - 1) + 1] = (y[0] + y[1]) / 2;

	  meml_matrix_element_set_f (restriction, mesh->edges[i][0] - 1,
				     zeiger - 1, 0.5);
	  meml_matrix_element_set_f (restriction, mesh->edges[i][1] - 1,
				     zeiger - 1, 0.5);

	  /* ist es bei randpunkt ? */
	  for (j = 0; j < mesh->number_of_boundary_edges; j++)
	    {
	      if (mesh->boundary_edges[j] == i)
		{
		  boundary->data[bzeiger] = zeiger;
		  bzeiger++;
		  break;
		}
	    }
	  zeiger++;
	}
    }

  /* how many new triangles */

  wieviele = 0;
  for (i = 0; i < number_of_triangles; i++)
    {
      /* wieviele kanten sind makiert */
      number_of_mark_edges = 0;
      for (j = 0; j < 3; j++)
	{
	  kante = mesh->triangle2edge[i][j];
	  if (mark_edges->data[kante] > 0)
	    number_of_mark_edges++;
	}

      /* if nothing had changes copy the old triangle */
      if (number_of_mark_edges == 0)
	{
	  /* nothing */
	  wieviele += 1;
	}

      /* If all three sides are divided, 
         new triangles are formed by joining the side midpoints. */

      if (number_of_mark_edges == 3)
	{
	  /* rot */
	  wieviele += 4;
	}

      /* If two sides are divided, 
         the midpoint of the longest edge is joined with the the 
         opposing corner and with the other midpoint. */

      if (number_of_mark_edges == 2)
	{
	  /* blau */
	  wieviele += 3;
	}

      /* If only the longest edge is divided, its midpoint is joined with the
         opposing corner. */

      if (number_of_mark_edges == 1)
	{
	  /* gruen */
	  wieviele += 2;
	}
    }

  triangles = meml_indexarray_new (3 * wieviele);

  /* Form the new triangles. */
  zeiger = 0;
  for (i = 0; i < number_of_triangles; i++)
    {

      number_of_mark_edges = 0;
      laengste = 0;
      welche = -1;

      for (j = 0; j < 3; j++)
	{
	  kante = mesh->triangle2edge[i][j];

	  kanten[j][0] = mesh->edges[kante][0];
	  kanten[j][1] = mesh->edges[kante][1];

	  /* wieviele kanten sind makiert ? */
	  /* moegliche mittelpunkte speichern 
	     die mittelpunkte muessen fuer die rote verfeinerung die gleiche 
	     nummerierung haben wie die kanten. 
	     mittelpunkte mit einer 0 sind keine echten, also ggf. abfangen */
	  if (mark_edges->data[kante] > 0)
	    {
	      mittelpunkte[j] = mark_edges->data[kante];
	      number_of_mark_edges++;
	    }
	  else
	    mittelpunkte[j] = 0;

	  /* dreieck zwischenspeichern */
	  t[j] = mesh->triangles->data[3 * i + j];

	  /* welche kante ist die laengste */
	  x[0] = mesh->points->data[2 * (mesh->edges[kante][0] - 1)];
	  x[1] = mesh->points->data[2 * (mesh->edges[kante][1] - 1)];
	  y[0] = mesh->points->data[2 * (mesh->edges[kante][0] - 1) + 1];
	  y[1] = mesh->points->data[2 * (mesh->edges[kante][1] - 1) + 1];

	  laenge =
	    (x[0] - x[1]) * (x[0] - x[1]) + (y[0] - y[1]) * (y[0] - y[1]);

	  if (laenge > laengste)
	    {
	      laengste = laenge;
	      welche = j;
	    }
	}

      /* welcher Punkt liegt der laengsten Kante gegenueber */
      for (j = 0; j < 3; j++)
	{
	  if ((t[j] != kanten[welche][0]) && (t[j] != kanten[welche][1]))
	    {
	      gegenueber = t[j];
	    }
	}

      /* if nothing had changes copy the old triangle */
      if (number_of_mark_edges == 0)
	{
	  /* nothing */
	  for (j = 0; j < 3; j++)
	    triangles->data[3 * zeiger + j] =
	      mesh->triangles->data[i * 3 + j];

	  mesh->son[i][0] = zeiger;
	  mesh->son[i][1] = -1;
	  mesh->son[i][2] = -1;
	  mesh->son[i][3] = -1;

	  zeiger += 1;
	}

      /* If all three sides are divided, 
         new triangles are formed by joining the side midpoints. */

      /* rot */
      else if (number_of_mark_edges == 3)
	{
	  /* erst das mittlere dreieck */
	  for (j = 0; j < 3; j++)
	      triangles->data[3 * zeiger + j] = mittelpunkte[j];
	  mesh->son[i][0] = zeiger;
	  zeiger += 1;

	  /* jetzt die, die auch alte punkte enthalten */
	  for (k = 0; k < 3; k++)
	    {
	      triangles->data[zeiger * 3] = t[k];
	      zaehler = 1;
	      for (j = 0; j < 3; j++)
		{
		  if ((t[k] == kanten[j][0]) || (t[k] == kanten[j][1]))
		    {
		      triangles->data[zeiger * 3 + zaehler] = mittelpunkte[j];
		      zaehler++;
		    }
		}
	      mesh->son[i][k+1] = zeiger;
	      zeiger += 1;
	    }
	}
      /* If two sides are divided, 
         the midpoint of the longest edge is joined with the the 
         opposing corner and with the other midpoint. */
      else if (number_of_mark_edges == 2)
	{
	  /* blau */
	  for (j = 0; j < 3; j++)
	    {
	      if ((mittelpunkte[j] > 0) && (welche != j))
		{
		  /* 1 dreieck */
		  triangles->data[zeiger * 3] = kanten[j][0];
		  triangles->data[zeiger * 3 + 1] = mittelpunkte[j];
		  triangles->data[zeiger * 3 + 2] = mittelpunkte[welche];
		  mesh->son[i][0] = zeiger;
		  zeiger += 1;

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

	      if (mittelpunkte[j] == 0)
		{
		  /* 3 dreieck */
		  triangles->data[zeiger * 3] = kanten[j][0];
		  triangles->data[zeiger * 3 + 1] = kanten[j][1];
		  triangles->data[zeiger * 3 + 2] = mittelpunkte[welche];
		  mesh->son[i][2] = zeiger;
		  zeiger += 1;
		}