qfel.c

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

  static const MEML_FLOAT S4[6][6] = {
    {2 / 24.0, 1 / 24.0, 1 / 24.0, 1.0 / 15, 1.0 / 30, 1.0 / 15},
    {1 / 24.0, 2 / 24.0, 1 / 24.0, 1.0 / 15, 1.0 / 15, 1.0 / 30},
    {1 / 24.0, 1 / 24.0, 2 / 24.0, 1.0 / 30, 1.0 / 15, 1.0 / 15},
    {1.0 / 15, 1.0 / 15, 1.0 / 30, 4.0 / 45, 2.0 / 45, 2.0 / 45},
    {1.0 / 30, 1.0 / 15, 1.0 / 15, 2.0 / 45, 4.0 / 45, 2.0 / 45},
    {1.0 / 15, 1.0 / 30, 1.0 / 15, 2.0 / 45, 2.0 / 45, 4.0 / 45}
  };

  MEML_INT i;
  MEML_INT x, y;
  MEML_FLOAT det;
  TRIANGLE the_triangle;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT *triangle = triangles->data;
  MEML_INT knoten[6];
  VECTOR * b;
  
  b = meml_vector_new(f->dim);

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

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

      /* 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++)
	knoten[x] = the_triangle.p[x] - 1;

      for (x = 0; x < 3; x++)
	knoten[x + 3] = mesh->triangles_adof[i][x] + mesh->number_of_points;

      for (x = 0; x < 6; x++)
	for (y = 0; y < 6; y++)
	  {
	      b->data[knoten[x]] += det * S4[x][y] * f->data[ knoten[y] ];
	  }
    }

  return (b);
}

/**
 * \~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 qfel_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);
      b->data[b_points[i] - 1] = value;
    }

  for (i = 0; i < mesh->number_of_boundary_edges; i++)
    {
      /* die boundary_edges entsprechen in ihrer nummerierung den 
         adof daher kann man einfach die boundary_edges nehmen um 
         heraus zu finden, welche der adof am rand liegen */
      meml_matrix_set_row_2_unit_vector (G,
					 mesh->number_of_points +
					 mesh->boundary_edges[i]);
      meml_matrix_set_col_2_unit_vector (G,
					 mesh->number_of_points +
					 mesh->boundary_edges[i]);
      b->data[mesh->number_of_points + mesh->boundary_edges[i]] = 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<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 qfel_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;
    VECTOR * hvalue;

    if (b !=NULL)
    {
	hvalue = qfel_nodal2hierarchical_base (value,mesh);

	for (i = 0; i < number_of_b_points; i++)
	{
	    if ( where->data[b_points[i] - 1] == 1 )
	    {
		b->data[b_points[i] - 1] = hvalue->data[b_points[i] - 1];
	    }
	}
	for (i = 0; i < mesh->number_of_boundary_edges; i++)
	{
	    /* die boundary_edges entsprechen in ihrer nummerierung den 
	       adof daher kann man einfach die boundary_edges nehmen um 
	       heraus zu finden, welche der adof am rand liegen */
	    if ( where->data[b_points[i] - 1] == 1 )
	    {
		b->data[mesh->number_of_points + mesh->boundary_edges[i]] = 
		    hvalue->data[mesh->number_of_points + mesh->boundary_edges[i]];
	    }
	}
	
	meml_vector_free(hvalue);
    }

    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);
	    }
	}
	for (i = 0; i < mesh->number_of_boundary_edges; i++)
	{
	    meml_matrix_set_row_2_unit_vector (G,
					       mesh->number_of_points +
					       mesh->boundary_edges[i]);
	}
    }
}

void qfel_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;

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

VECTOR *qfel_hierarchical2nodal_base (const VECTOR * x, const MESH * mesh)
{
  VECTOR *x_copy;

  x_copy = meml_vector_copy (x);

  qfel_hierarchical2nodal_base_f (x_copy, mesh);

  return (x_copy);
}


void qfel_hierarchical2nodal_base_f (VECTOR * x, const MESH * mesh)
{
  ME_MATRIX *T;

  T = qfel_hierarchical2nodal_base_ff (mesh);

  meml_matrix_vector_mul_ff (T, x);

  meml_matrix_free (T);
}

ME_MATRIX *qfel_hierarchical2nodal_base_ff (const MESH * mesh)
{
  ME_MATRIX *B, *R;
  const MEML_INT size = mesh->number_of_points + mesh->number_of_adof;
  MEML_INT i, x, y;


  B = meml_matrix_new_eye (LS, size, size);

  for (i = 0; i < mesh->number_of_edges; i++)
    {
      x = mesh->edges[i][0] - 1;
      y = mesh->edges[i][1] - 1;
      meml_matrix_element_set (B, i + mesh->number_of_points, x, 0.5);
      meml_matrix_element_set (B, i + mesh->number_of_points, y, 0.5);
    }

  R = meml_matrix_convert (CS, B);
  meml_matrix_free (B);

  return (R);
}


VECTOR *qfel_nodal2hierarchical_base (const VECTOR * x, const MESH * mesh)
{
  VECTOR *x_copy;

  x_copy = meml_vector_copy (x);

  qfel_nodal2hierarchical_base_f (x_copy, mesh);

  return (x_copy);
}


void qfel_nodal2hierarchical_base_f (VECTOR * x, const MESH * mesh)
{
  ME_MATRIX *T;

  T = qfel_nodal2hierarchical_base_ff (mesh);

  meml_matrix_vector_mul_ff (T, x);

  meml_matrix_free(T);
}

ME_MATRIX *qfel_nodal2hierarchical_base_ff (const MESH * mesh)
{
  ME_MATRIX *B, *R;
  const MEML_INT size = mesh->number_of_points + mesh->number_of_adof;
  MEML_INT i, x, y;

  B = meml_matrix_new_eye (LS, size, size);

  for (i = 0; i < mesh->number_of_edges; i++)
    {
      x = mesh->edges[i][0] - 1;
      y = mesh->edges[i][1] - 1;
      meml_matrix_element_set (B, i + mesh->number_of_points, x, -0.5);
      meml_matrix_element_set (B, i + mesh->number_of_points, y, -0.5);
    }

  R = meml_matrix_convert (CS, B);
  meml_matrix_free (B);

  return (R);
}


/**
 * \~german 
 * \brief Assembliert die rechte Seite des LGS fuer das Galerkinverfahren
 * 
 * Der Vector f wird in der nodalen Form erwartet.
 */
double qfel_integrating (const VECTOR * f, MESH * mesh)
{
  INDEXARRAY *triangles = mesh->triangles;
  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[6][7] = {
    {Q1 (gaus_x[0], gaus_y[0]), Q1 (gaus_x[1], gaus_y[1]),
     Q1 (gaus_x[2], gaus_y[2]), Q1 (gaus_x[3], gaus_y[3]),
     Q1 (gaus_x[4], gaus_y[4]), Q1 (gaus_x[5], gaus_y[5]),
     Q1 (gaus_x[6], gaus_y[6])},
    {Q2 (gaus_x[0], gaus_y[0]), Q2 (gaus_x[1], gaus_y[1]),
     Q2 (gaus_x[2], gaus_y[2]), Q2 (gaus_x[3], gaus_y[3]),
     Q2 (gaus_x[4], gaus_y[4]), Q2 (gaus_x[5], gaus_y[5]),
     Q2 (gaus_x[6], gaus_y[6])},
    {Q3 (gaus_x[0], gaus_y[0]), Q3 (gaus_x[1], gaus_y[1]),
     Q3 (gaus_x[2], gaus_y[2]), Q3 (gaus_x[3], gaus_y[3]),
     Q3 (gaus_x[4], gaus_y[4]), Q3 (gaus_x[5], gaus_y[5]),
     Q3 (gaus_x[6], gaus_y[6])},
    {Q4 (gaus_x[0], gaus_y[0]), Q4 (gaus_x[1], gaus_y[1]),
     Q4 (gaus_x[2], gaus_y[2]), Q4 (gaus_x[3], gaus_y[3]),
     Q4 (gaus_x[4], gaus_y[4]), Q4 (gaus_x[5], gaus_y[5]),
     Q4 (gaus_x[6], gaus_y[6])},
    {Q5 (gaus_x[0], gaus_y[0]), Q5 (gaus_x[1], gaus_y[1]),
     Q5 (gaus_x[2], gaus_y[2]), Q5 (gaus_x[3], gaus_y[3]),
     Q5 (gaus_x[4], gaus_y[4]), Q5 (gaus_x[5], gaus_y[5]),
     Q5 (gaus_x[6], gaus_y[6])},
    {Q6 (gaus_x[0], gaus_y[0]), Q6 (gaus_x[1], gaus_y[1]),
     Q6 (gaus_x[2], gaus_y[2]), Q6 (gaus_x[3], gaus_y[3]),
     Q6 (gaus_x[4], gaus_y[4]), Q6 (gaus_x[5], gaus_y[5]),
     Q6 (gaus_x[6], gaus_y[6])}
  };

  MEML_INT i, x, y;
  MEML_FLOAT det;
  TRIANGLE the_triangle;
  MEML_FLOAT I_temp[6];
  MEML_FLOAT Integral = 0;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT *triangle = triangles->data;
  MEML_FLOAT f_value[6];
  MEML_INT knoten[6];
  VECTOR *f_copy;

  f_copy = qfel_nodal2hierarchical_base (f, mesh);

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

      for (x = 0; x < 3; x++)
	{
	  knoten[x] = the_triangle.p[x] - 1;
	}

      for (x = 0; x < 3; x++)
	{
	  knoten[x + 3] = mesh->triangles_adof[i][x] + mesh->number_of_points;
	}

      for (x = 0; x < 6; x++)
	{
	  f_value[x] = f_copy->data[knoten[x]];
	}

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

      for (x = 0; x < 6; x++)
	{
	  I_temp[x] = 0;
	  for (y = 0; y < 7; y++)
	    {
	      I_temp[x] += phi[x][y] * f_value[x] * gaus_w[y];
	    }
	  I_temp[x] = I_temp[x] * det;
	}

      for (x = 0; x < 6; x++)
	{
	  Integral += I_temp[x];
	}
    }

  return (Integral);
}


/**
 * \brief applies the values of the function f to the vector x
 *
 * The order in x depends on the order of the points in mesh. 
 * t is the time argument of the function f. 
 *  
 */
void qfel_apply_function_to_vector_f (VECTOR * v, MESH * mesh,
				      xyt_function f, MEML_FLOAT t)
{
  MEML_INT i;
  MEML_FLOAT x, y;

  for (i = 0; i < mesh->number_of_points; i++)
    {
      x = mesh->points->data[2 * i];
      y = mesh->points->data[(2 * i) + 1];

      v->data[i] = (*f) (x, y, t);
    }

  for (i = 0; i < mesh->number_of_adof; i++)
    {
      x = mesh->additional_degrees_of_freedom->data[2 * i];
      y = mesh->additional_degrees_of_freedom->data[(2 * i) + 1];

      v->data[i + mesh->number_of_points] = (*f) (x, y, t);
    }
}

/**
 * \brief applies the values of the function f to the vector x
 *
 * The order in x depends on the order of the points in mesh. 
 * t is the time argument of the function f. 
 *  
 */
VECTOR *qfel_apply_function_to_vector (MESH * mesh, xyt_function f,
				       MEML_FLOAT t)
{
  VECTOR *v;
  MEML_INT size = mesh->number_of_points + mesh->number_of_adof;

  v = vecl_vector_new (size);

  qfel_apply_function_to_vector_f (v, mesh, f, t);

  return (v);
}


/**
 * \brief 
 * \~german 
 *  Berechnet den Gradienten einer Funktion, die durch ihre Werte an 
 *  den Knoten stellen des Gitters gegeben ist. 
 * \~english
 *  Computes the gradient of a function, which is given by the 
 *  values of the function at the nodes of the mesh.
 *
 * \f$ (v_x , v_y) = grad \; v  \f$<br>
 *
 * \~german