lfel.c

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

 */
int lfel_assemb_diffusion_f (const VECTOR * eps, ME_MATRIX * A, MESH * mesh)
{
  INDEXARRAY *triangles = mesh->triangles;
  VECTOR *points = mesh->points;
  TRANSFORMATION *transformation = mesh->transformation;

  const MEML_FLOAT S1[3][3] = {
    {0.5, -0.5, 0}, {-0.5, 0.5, 0}, {0, 0, 0}
  };
  const MEML_FLOAT S2[3][3] = {
    {1, -0.5, -0.5}, {-0.5, 0, 0.5}, {-0.5, 0.5, 0}
  };
  const MEML_FLOAT S3[3][3] = {
    {0.5, 0, -0.5}, {0, 0, 0}, {-0.5, 0, 0.5}
  };

  MEML_INT i, row, col;
  int x, y;
  MEML_FLOAT det, gamma1, gamma2, gamma3, temp, epsontriangle;
  TRIANGLE the_triangle;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT *triangle = triangles->data;
  M_TYPE type;


  if ( (eps->dim != 1) && (eps->dim == number_of_triangles) 
       && (eps->dim != mesh->number_of_triangles))
    {
      fprintf (stderr, "lfel_assemb_diffusion : "
	       "The size of the eps indexarray does not fit!\n");
      return (EXIT_FAILURE);
    }
    
  meml_matrix_property_get (A, &col, &row, &type);

  if ((row < points->dim / 2) || (col < points->dim / 2))
    {
      fprintf (stderr, "lfel_assemb_diffusion : "
	       "The size of the matrix is not big enough "
	       "so store the data!\n");
      return (EXIT_FAILURE);
    }

  for (i = 0; i < number_of_triangles; i++)
    {
      the_triangle = smfl_get_triangle (i, triangle);

      det = fabs (transformation[i].det);

      gamma1 = (transformation[i].B[0][1] * transformation[i].B[0][1]
		+
		transformation[i].B[1][1] * transformation[i].B[1][1]) / det;

      gamma2 = -(transformation[i].B[0][0] * transformation[i].B[0][1]
		 +
		 transformation[i].B[1][0] * transformation[i].B[1][1]) / det;

      gamma3 = (transformation[i].B[1][0] * transformation[i].B[1][0] +
		transformation[i].B[0][0] * transformation[i].B[0][0]) / det;

      /* Warnung. MatLab nutzt 1,1 als obersten Eintrag
         C hingegen 0,0. Das ist ein gut moeglichkeit
         speicherzugriffsfehler zu produzieren */

      if (eps->dim == 1)
	{
	  epsontriangle = eps->data[0];
	}
      else if (eps->dim == number_of_triangles)
	{
	  epsontriangle = eps->data[i];
	}
      else 
	{
	  epsontriangle = 0;
	  for (x=0;x<3;x++)
	    {
	      epsontriangle += eps->data[the_triangle.p[x] - 1];
	    }
	  epsontriangle = epsontriangle /3;
	}

      for (x = 0; x < 3; x++)
	for (y = 0; y < 3; y++)
	  {
	    temp =
	      epsontriangle * (gamma1 * S1[x][y] + gamma2 * S2[x][y] +
		     gamma3 * S3[x][y]);

	    /* wegen C <-> Matlab -1 */
	    meml_matrix_element_add_f (A, the_triangle.p[x] - 1,
				       the_triangle.p[y] - 1, temp);
	  }
    }

  return (EXIT_SUCCESS);
}

/**
 * \~german 
 *  \brief Setzt homogene Dirichlet-Randbedingungen 
 *
 *   Sie Funktion ersetzt Zeilen und Spalten fuer die Dirichelt-Randbedingungen 
 *   gelten durch die entsprechenden Einheitsvektoren. Auf der rechten Seite 
 *   werden die entsprechenden Stellen durch value ersetzt<br>
 *   <br><b>
 *   Mittelfristig sollte eine weitere Funktion geschrieben werden, 
 *   die stattdessen ein reduziertes Gleichungssystem berechnet.</b><br>
 *  
 *  \param G fertig assemblierte Galerkin-Matrix
 *  \param b die rechte Seite des LGS
 *  \param value der gesuchten Funktion am Rand
 */
void lfel_set_h_dirichlet_boundary_cond (ME_MATRIX * G, VECTOR * b,
					 MEML_FLOAT value, MESH * mesh)
{
  INDEXARRAY *boundary = mesh->boundary;
  MEML_INT i;
  MEML_INT number_of_b_points = boundary->dim;
  MEML_INT *b_points = boundary->data;

  for (i = 0; i < number_of_b_points; i++)
    {
      meml_matrix_set_row_2_unit_vector (G, b_points[i] - 1);
      meml_matrix_set_col_2_unit_vector (G, b_points[i] - 1);
    }
  
  if (b!=NULL)
    for (i = 0; i < number_of_b_points; i++)
      {
	b->data[b_points[i] - 1] = value;
      }
}

/**
 * \~german 
 *  \brief Setzt homogene Dirichlet-Randbedingungen 
 *
 *   Sie Funktion ersetzt Zeilen und Spalten fuer die Dirichelt-Randbedingungen 
 *   gelten durch die entsprechenden Einheitsvektoren. Auf der rechten Seite 
 *   werden die entsprechenden Stellen durch value ersetzt, falls der 
 *   entsprechende indexarray eintrag auf 1 steht.<br>
 */
void lfel_set_dirichlet_boundary_cond (ME_MATRIX * G, VECTOR * b,VECTOR * value,
				       MESH * mesh)
{
  INDEXARRAY * where;
  MEML_INT i;

  where=meml_indexarray_new(value->dim);

  for (i=0;i<where->dim;i++)
    where->data[i]=1;

  lfel_set_dirichlet_boundary_cond_f (G,b,where,value,mesh);
  
  meml_indexarray_free(where);
}


/**
 * \~german 
 *  \brief Setzt homogene Dirichlet-Randbedingungen 
 *
 *   Sie Funktion ersetzt Zeilen und Spalten fuer die Dirichelt-Randbedingungen 
 *   gelten durch die entsprechenden Einheitsvektoren. Auf der rechten Seite 
 *   werden die entsprechenden Stellen durch value ersetzt, falls der 
 *   entsprechende indexarray eintrag auf 1 steht.<br>
 */
void lfel_set_dirichlet_boundary_cond_f (ME_MATRIX * G, VECTOR * b,
					 INDEXARRAY * where, VECTOR * value,
					 MESH * mesh)
{
  INDEXARRAY *boundary = mesh->boundary;
  MEML_INT i;
  MEML_INT number_of_b_points = boundary->dim;
  MEML_INT *b_points = boundary->data;

  if (b !=NULL)
    for (i = 0; i < number_of_b_points; i++)
      {
	if ( where->data[b_points[i] - 1] == 1 )
	  {
	    b->data[b_points[i] - 1] = value->data[b_points[i] - 1];
	  }
      }
  if (G !=NULL)
    for (i = 0; i < number_of_b_points; i++)
      {
	if ( where->data[b_points[i] - 1] == 1 )
	  {
	    meml_matrix_set_row_2_unit_vector (G, b_points[i] - 1);
	  }
      }
}

INDEXARRAY * lfel_return_nodes(MESH * mesh, xyt_function f, MEML_FLOAT t)
{
  MEML_INT i;
  MEML_FLOAT x, y, z;
  INDEXARRAY * which;
  
  which = meml_indexarray_new(mesh->number_of_points);

  for (i = 0; i < mesh->number_of_points; i++)
    {
      x = mesh->points->data[2 * i];
      y = mesh->points->data[(2 * i) + 1];
      
      z = f(x,y,t);
      
      if ( fabs(z) < MEML_EPS )
      {
	which->data[i] = 1;
      }
    }
  
  return(which);
}

int lfel_assemb_reaction (const MEML_FLOAT r, ME_MATRIX * A, MESH * mesh)
{
  VECTOR * rarray;
  int re;

  rarray = meml_vector_new(1);

  rarray->data[0]=r;

  re=lfel_assemb_reaction_f (rarray, A, mesh);

  meml_vector_free(rarray);
  
  return(re);
}

/**
 * \~german 
 *  \brief Assembliert die Massen-Matrix/den Reaktionsanteil
 *   fuer das Galerkinverfahren
 *
 *   Die Matrix A muss vor dem Aufruf dieser Funktion erzeugt worden sein.
 *   Bei groesseren Gitter sollte A eine Sparse-Matrix-Typ sein.
 *  
 *  \param r Reaktionsparameter der partiellen Differentialgleichung
 *  \param A die Matrix zu der die Steifigkeitsmatrix addiert werden soll
 */
int lfel_assemb_reaction_f (const VECTOR * r, ME_MATRIX * A, MESH * mesh)
{
  INDEXARRAY *triangles = mesh->triangles;
  VECTOR *points = mesh->points;
  TRANSFORMATION *transformation = mesh->transformation;

  static const MEML_FLOAT S4[3][3] = {
    {2, 1, 1}, {1, 2, 1}, {1, 1, 2}
  };

  MEML_INT i, row, col;
  MEML_INT x, y;
  MEML_FLOAT det, gamma, temp;
  TRIANGLE the_triangle;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT *triangle = triangles->data;
  M_TYPE type;


  if ( (r->dim != 1) && (r->dim == number_of_triangles) 
       && (r->dim != mesh->number_of_triangles))
    {
      fprintf (stderr, "lfel_assemb_reaction : "
	       "The size of the eps indexarray does not fit!\n");
      return (EXIT_FAILURE);
    }

  meml_matrix_property_get (A, &col, &row, &type);

  if ((row < points->dim / 2) || (col < points->dim / 2))
    {
      fprintf (stderr, "The size of the matrix is not big enough "
	       "so store the data!\n");
      return (EXIT_FAILURE);
    }

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

      the_triangle = smfl_get_triangle (i, triangle);

      det = fabs (transformation[i].det);

      if (r->dim == 1)
	{
	  gamma = (det / 24) * r->data[0];
	}
      else if (r->dim == number_of_triangles)
	{
	  gamma = (det / 24) * r->data[i];
	}
      else 
	{
	  temp = 0;
	  for (x=0;x<3;x++)
	    temp += r->data[the_triangle.p[x] - 1];
	  temp = temp /3;
	  
	  gamma = (det / 24) * temp;
	}

      /* Warnung. MatLab nutzt 1,1 als obersten Eintrag
         C hingegen 0,0. Das ist ein gut moeglichkeit
         speicherzugriffsfehler zu produzieren */


      for (x = 0; x < 3; x++)
	for (y = 0; y < 3; y++)
	  {
	    temp = gamma * S4[x][y];
	    meml_matrix_element_add (A, the_triangle.p[x] - 1,
				     the_triangle.p[y] - 1, temp);
	  }
    }

  return (EXIT_SUCCESS);
}

/**
 * \brief
 * \~english calculates the the right hand side  
 * \~german Assembliert die rechte Seite des LGS fuer das Galerkinverfahren
 *
 *  \param f Vektor der die Werte der Funktion f an jedem Knoten enthaelt. 
 *   Die Reihenfolge der Daten in f muss der Nummerierung der Knoten 
 *   entsprechen.
 *   Die Dimension von f entspricht folglich der Anzahl der Knoten.
 */
VECTOR *lfel_assemb_right_side (const VECTOR * f, const MESH * mesh)
{
  VECTOR *r, *temp;

  r = meml_vector_new_ones(mesh->number_of_triangles);
  temp = lfel_assemb_right_side_f (f, r, mesh);
  meml_vector_free(r);

  return(temp);
}


VECTOR *lfel_assemb_right_side_f (const VECTOR * f, const VECTOR *r, const MESH * mesh)
{
  INDEXARRAY *triangles = mesh->triangles;
  VECTOR *points = mesh->points;
  TRANSFORMATION *transformation = mesh->transformation;

  const MEML_FLOAT gaus_x[3] = { 0.5, 0.0, 0.5 };
  const MEML_FLOAT gaus_y[3] = { 0.0, 0.5, 0.5 };
  const MEML_FLOAT gaus_w[3] = { 1.0 / 6.0, 1.0 / 6.0, 1.0 / 6.0 };

  const MEML_FLOAT phi[3][3] = {
    {N1 (gaus_x[0], gaus_y[0]), N1 (gaus_x[1], gaus_y[1]),
     N1 (gaus_x[2], gaus_y[2])},
    {N2 (gaus_x[0], gaus_y[0]), N2 (gaus_x[1], gaus_y[1]),
     N2 (gaus_x[2], gaus_y[2])},
    {N3 (gaus_x[0], gaus_y[0]), N3 (gaus_x[1], gaus_y[1]),
     N3 (gaus_x[2], gaus_y[2])}
  };


  MEML_INT i, x;
  MEML_FLOAT temp, det;
  TRIANGLE the_triangle;
  MEML_FLOAT f_temp[3];
  MEML_INT number_of_points = points->dim / 2;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT *triangle = triangles->data;
  MEML_FLOAT *b;
  VECTOR *temp_vector;
  MEML_FLOAT f_value[3];
  MEML_FLOAT rOnTriangle;

  if (f->dim != number_of_points)
    {
      fprintf (stderr, "lfel_assemb_right_side : "
	       "Dimension of f does not agree with the number of points");
      exit (EXIT_FAILURE);
    }

  temp_vector = vecl_vector_new (number_of_points);

  b = temp_vector->data;

  for (i = 0; i < number_of_points; i++)
    b[i] = 0;

  for (i = 0; i < number_of_triangles; i++)
    {
      the_triangle = smfl_get_triangle (i, triangle);


      if (r->dim == 1)
	{
	  rOnTriangle = r->data[0];
	}
      else if (r->dim == number_of_triangles)
	{
	  rOnTriangle = r->data[i];
	}
      else 
	{
	  temp = 0;
	  for (x=0;x<3;x++)
	    temp += r->data[the_triangle.p[x] - 1];
	  temp = temp /3;

	  rOnTriangle = temp;
	}

      /* den Wert von f an den Mittelpunkten der Seiten zu berechnen */
      f_value[0] =
	(f->data[the_triangle.p[0] - 1] + f->data[the_triangle.p[1] - 1]) / 2;
      f_value[1] =
	(f->data[the_triangle.p[0] - 1] + f->data[the_triangle.p[2] - 1]) / 2;
      f_value[2] =
	(f->data[the_triangle.p[1] - 1] + f->data[the_triangle.p[2] - 1]) / 2;

      f_temp[0] = 0;
      f_temp[1] = 0;
      f_temp[2] = 0;

      for (x = 0; x < 3; x++)
	{
	  temp = gaus_w[x] * f_value[x] * rOnTriangle;
	  f_temp[0] = f_temp[0] + temp * phi[0][x];
	  f_temp[1] = f_temp[1] + temp * phi[1][x];
	  f_temp[2] = f_temp[2] + temp * phi[2][x];
	}

      det = fabs (transformation[i].det);

      for (x = 0; x < 3; x++)
	{
	  b[the_triangle.p[x] - 1] =
	    b[the_triangle.p[x] - 1] + f_temp[x] * det;
	}
    }

  return (temp_vector);
}

/**
 * \~german 
 *  \brief Assembliert die rechte Seite des LGS fuer das Galerkinverfahren
 */
VECTOR *lfel_assemb_load_vector (xyt_function f, const double time,
				 const MESH * mesh)
{
    return(
	lfel_assemb_load_vector_f (f, time, NULL, mesh) );
}


VECTOR *lfel_assemb_load_vector_f (xyt_function f, const double time, VECTOR * r, 
				   const MESH * mesh)
{
  INDEXARRAY *triangles = mesh->triangles;
  VECTOR *points = mesh->points;
  TRANSFORMATION *transformation = mesh->transformation;

  const MEML_FLOAT gaus_x[7] = { 0, 1, 0, 0.5, 0.5, 0, 1.0 / 3.0 };
  const MEML_FLOAT gaus_y[7] = { 0, 0, 1, 0, 0.5, 0.5, 1.0 / 3.0 };
  const MEML_FLOAT gaus_w[7] = { 1.0 / 40.0, 1.0 / 40.0, 1.0 / 40.0,
    1.0 / 15.0, 1.0 / 15.0, 1.0 / 15.0,
    27.0 / 120.0
  };

  const MEML_FLOAT phi[3][7] = {
    {N1 (gaus_x[0], gaus_y[0]), N1 (gaus_x[1], gaus_y[1]),
     N1 (gaus_x[2], gaus_y[2]), N1 (gaus_x[3], gaus_y[3]),
     N1 (gaus_x[4], gaus_y[4]), N1 (gaus_x[5], gaus_y[5]),
     N1 (gaus_x[6], gaus_y[6])},
    {N2 (gaus_x[0], gaus_y[0]), N2 (gaus_x[1], gaus_y[1]),
     N2 (gaus_x[2], gaus_y[2]), N2 (gaus_x[3], gaus_y[3]),
     N2 (gaus_x[4], gaus_y[4]), N2 (gaus_x[5], gaus_y[5]),
     N2 (gaus_x[6], gaus_y[6])},
    {N3 (gaus_x[0], gaus_y[0]), N3 (gaus_x[1], gaus_y[1]),
     N3 (gaus_x[2], gaus_y[2]), N3 (gaus_x[3], gaus_y[3]),
     N3 (gaus_x[4], gaus_y[4]), N3 (gaus_x[5], gaus_y[5]),
     N3 (gaus_x[6], gaus_y[6])}
  };

  MEML_INT i, x;
  MEML_FLOAT temp, det;
  TRIANGLE the_triangle;
  MEML_FLOAT f_temp[3];
  MEML_INT number_of_points = points->dim / 2;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT *triangle = triangles->data;
  MEML_FLOAT *b;
  VECTOR *temp_vector;
  MEML_FLOAT f_value[7];
  MEML_FLOAT d[2];
  VECTOR * triangleR;
  temp_vector = meml_vector_new (number_of_points);

  b = temp_vector->data;

  for (i = 0; i < number_of_points; i++)
    b[i] = 0;


  if (r != NULL)
  {
      if ( r->dim == mesh->number_of_triangles)
	  triangleR = r; 
      else if ( r->dim == 1)
      {
	  triangleR = meml_vector_new_ones(mesh->number_of_triangles);      
	  meml_vector_scaling_f