qfel.c

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

 * \param v sind die Werte einer Funktion an den Knoten der 
 * Triangulierung. Sie werden in einem Vektor in der Reihenfolge der 
 * Knoten abgespeichert.
 *
 *  \param transformation Abildungen des Referenzdreiecks auf die Dreiecke der 
 *   Triangulierung lfel_compute_transformations
 *
 *  \param triangles enthaelt die Informationen ueber die Knoten der Dreiecke 
 *  der Triangulierung in einem Format wie es von der Funktion 
 *  lfel_read_triangle_data erzeugt wird.
 *
 *  \param v_x enthaelt nach der Beendigung der Funktion
 *  die partielle Ableitung von v in x Richtung
 *
 *  \param v_y enthaelt nach der Beendigung der Funktion
 *  die partielle Ableitung von v in y Richtung
 *
 *  v_x und v_y muessen vor dem Aufruf der Funktion initialisiert worden sein.
 *
 */
void qfel_grad (const VECTOR * v, MESH * mesh, VECTOR * v_x, VECTOR * v_y)
{
  TRANSFORMATION *transformation = mesh->transformation;
  MEML_INT number_of_triangles = mesh->number_of_triangles;

  MEML_FLOAT I[2][2];
  MEML_FLOAT lgrad[6][2];
  MEML_FLOAT temp[2], det, Omega;
  MEML_INT knoten[6];
  MEML_INT t, x, y, i, triangle, aktueller_knoten, vorzeichen;

  TRIANGLE the_triangle;
  VECTOR *vh;

  MEML_FLOAT tref_grad[6][6][2] = {
    /* 
       Punkte im Dreieick : 
       (0,0), (0,1), (1,0), (0.5 , 0) , (0.5, 0.5) , (0,0.5) 
     */
    /* grad Q1 = ( -1 , -1) */
    {
     {-1, -1}, {-1, -1}, {-1, -1}, {-1, -1}, {-1, -1}, {-1, -1}
     },
    /* grad Q2 = ( 1 , 0) */
    {
     {1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 0}, {1, 0}
     },
    /* grad Q3 = ( 0 , 1) */
    {
     {0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}, {0, 1}
     },
    /* grad Q4 = ( 4-8x-4y , -4x) */
    {
     {4, 0}, {0, 0}, {-4, -4}, {0, -2}, {-2, -2}, {2, 0}
     },
    /* grad Q5 = ( 4y , 4x) */
    {
     {0, 0}, {4, 0}, {0, 4}, {0, 2}, {2, 2}, {2, 0}
     },
    /* grad Q6 = ( -4y , 4-8y-4x ) */
    {
     {0, 4}, {-4, -4}, {0, 0}, {0, 2}, {-2, -2}, {-2, 0}
     }
  };

  vh = qfel_nodal2hierarchical_base (v, mesh);

  for (i = 0; i < v_x->dim; i++)
    {
      v_x->data[i] = 0;
      v_y->data[i] = 0;
    }

  for (t = 0; t < number_of_triangles; t++)
    {
      the_triangle = smfl_get_triangle (t, mesh->triangles->data);

      det = transformation[t].det;

      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[t][x] + mesh->number_of_points;

      for (x = 0; x < 6; x++)
	{
	  lgrad[x][0] = 0;
	  lgrad[x][1] = 0;
	}

      qfel_invtrans_22_Matrix (mesh->transformation[t].B, I);

      for (x = 0; x < 6; x++)	/* Funktionen */
	{
	  for (y = 0; y < 6; y++)	/* Punkte     */
	    {
	      qfel_22mul (I, tref_grad[x][y], temp);
	      vorzeichen = (fabs (det) / det);
	      lgrad[y][0] += temp[0] * vh->data[knoten[x]] * vorzeichen;
	      lgrad[y][1] += temp[1] * vh->data[knoten[x]] * vorzeichen;
	    }
	}

      for (x = 0; x < 3; x++)
	{
	  /* kein /det wegen der inv-funktion */
	  v_x->data[knoten[x]] += lgrad[x][0] / 2;
	  v_y->data[knoten[x]] += lgrad[x][1] / 2;
	}
      for (x = 3; x < 6; x++)
	{
	  /* hier schon */
	  v_x->data[knoten[x]] += lgrad[x][0] / (2 * det) * vorzeichen;
	  v_y->data[knoten[x]] += lgrad[x][1] / (2 * det) * vorzeichen;
	}
    }

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

      v_x->data[aktueller_knoten] = v_x->data[aktueller_knoten] * 2;
      v_y->data[aktueller_knoten] = v_y->data[aktueller_knoten] * 2;
    }

  for (i = 0; i < mesh->number_of_points; i++)
    {
      Omega = 0;

      /* bereche die flaeche aller dreiecke die p als punkt haben */
      for (x = 0; x < mesh->point2triangle[i]->dim; x++)
	{
	  triangle = mesh->point2triangle[i]->data[x];

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

	  Omega += det / 2;
	}

      v_x->data[i] = v_x->data[i] / Omega;
      v_y->data[i] = v_y->data[i] / Omega;
    }

  meml_vector_free (vh);
}

#ifndef DOXYGEN_SHOULD_SKIP_THIS
inline void q_gibt_punkt (const VECTOR * u1, MESH * mesh,
			  MEML_INT node, MEML_FLOAT * x, MEML_FLOAT * y,
			  MEML_FLOAT * z)
{
  (*z) = u1->data[node];

  if (node < mesh->number_of_points)
    {
      (*x) = mesh->points->data[2 * (node)];
      (*y) = mesh->points->data[2 * (node) + 1];
    }
  else
    {
      (*x) =
	mesh->additional_degrees_of_freedom->data[2 *
						  (node -
						   mesh->number_of_points)];
      (*y) =
	mesh->additional_degrees_of_freedom->data[2 *
						  (node -
						   mesh->
						   number_of_points) + 1];
    }
}

INDEXARRAY *q_hole_indizes (MESH * mesh, MEML_INT node)
{
  MEML_INT anzahl_der_dreiecke;
  INDEXARRAY *liste, *welche_dreiecke, *liste2;
  TRIANGLE the_triangle;
  MEML_INT *triangle_node = mesh->triangles->data;
  MEML_INT i, j, zaehler = 0;


  /* wieviele dreiecke werden einbezogen */
  if (node < mesh->number_of_points)
    anzahl_der_dreiecke = mesh->point2triangle[node]->dim;
  else
    {
      if (mesh->edge2triangle[node - mesh->number_of_points][1] != -1)
	anzahl_der_dreiecke = 2;
      else
	anzahl_der_dreiecke = 1;
    }

  liste = meml_indexarray_new (anzahl_der_dreiecke * 6);
  welche_dreiecke = meml_indexarray_new (anzahl_der_dreiecke);

  /* welche dreiecke werden einbezogen */
  if (node < mesh->number_of_points)
    {
      for (i = 0; i < anzahl_der_dreiecke; i++)
	welche_dreiecke->data[i] = mesh->point2triangle[node]->data[i];
    }
  else
    {
      welche_dreiecke->data[0] =
	mesh->edge2triangle[node - mesh->number_of_points][0];
      if (anzahl_der_dreiecke == 2)
	welche_dreiecke->data[1] =
	  mesh->edge2triangle[node - mesh->number_of_points][1];
    }


  /* punkte aus den dreiecken holen */
  for (i = 0; i < anzahl_der_dreiecke; i++)
    {
      the_triangle =
	smfl_get_triangle (welche_dreiecke->data[i], triangle_node);

      /* erst die normalen punkte */
      for (j = 0; j < 3; j++)
	{
	  liste->data[zaehler] = the_triangle.p[j] - 1;

	  zaehler++;
	}
      /* dann die adof */
      for (j = 0; j < 3; j++)
	{
	  liste->data[zaehler] = mesh->number_of_points +
	    mesh->triangles_adof[welche_dreiecke->data[i]][j];

	  zaehler++;
	}
    }

  liste2 = doppelte_aussortieren (liste);

  meml_indexarray_free (liste);

  zaehler = 0;

  liste = meml_indexarray_new (liste2->dim - 1);

  /* den entwicklungspunkt ausfiltern */
  for (i = 0; i < liste2->dim; i++)
    {
      if (liste2->data[i] != (node))
	{
	  liste->data[zaehler] = liste2->data[i];
	  zaehler++;
	}
    }

  meml_indexarray_free (liste2);
  meml_indexarray_free (welche_dreiecke);

  return (liste);
}


VECTOR *apply_taylor (const MEML_INT node, INDEXARRAY * liste,
		      MESH * mesh, const VECTOR * u1)
{
  MEML_FLOAT entwx, entwy, entwz;
  MEML_FLOAT h1, h2, px, py, pz;
  VECTOR *b, *result;
  ME_MATRIX *A;
  MEML_INT i;

  /* daten des entwicklungspunktes holen */

  q_gibt_punkt (u1, mesh, node, &entwx, &entwy, &entwz);

  b = meml_vector_new (liste->dim);

  if (liste->dim < 12)
    {
      A = meml_matrix_new (ND, liste->dim, 5);
      result = meml_vector_new (5);
    }
  else
    {
      A = meml_matrix_new (ND, liste->dim, 9);
      result = meml_vector_new (9);
    }

  for (i = 0; i < liste->dim; i++)
    {
      q_gibt_punkt (u1, mesh, liste->data[i], &px, &py, &pz);

      h1 = px - entwx;
      h2 = py - entwy;

      b->data[i] = pz - entwz;

      /* f(a+h1,b+h2) = 
         f(a,b) + f_x h1 + f_y h2 + 
         0.5 f_xx h1^2 + 0.5 f_yy h2^2 + f_xy h1 h2 */
      sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 0, h1);
      sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 1, h2);
      sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 2, h1 * h1 / 2);
      sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 3, h2 * h2 / 2);
      sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 4, h1 * h2);
    }

  if (liste->dim > 11)
    {
      for (i = 0; i < liste->dim; i++)
	{
	  q_gibt_punkt (u1, mesh, liste->data[i], &px, &py, &pz);

	  h1 = px - entwx;
	  h2 = py - entwy;

	  b->data[i] = pz - entwz;

	  /* 
	     f(a+h1,b+h2) = 
	     f(a,b) + f_x h1 + f_y h2 + 
	     0.5 f_xx h1^2 + 0.5 f_yy h2^2 + f_xy h1 h2 
	     f_xxx h1^3 /6 + f_xxy h1^2*h2 /2 + 
	     f_xyy h1*h2^2 /2 + f_yyy h2^3 /6
	   */
	  sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 5,
				 h1 * h1 * h1 / 6);
	  sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 6,
				 h1 * h1 * h2 / 2);
	  sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 7,
				 h1 * h2 * h2 / 2);
	  sdml_nd_element_set_f (A->data.fortran_dense_matrix, i, 8,
				 h2 * h2 * h2 / 6);
	}
    }

  if (MEML_FALSE == itsl_least_square_solve_f (A, b, result))
    {
      printf ("liste dim %d \n", liste->dim);
      printf ("Panik bei punkt %d \n", node);
      meml_matrix_print(A);
      meml_vector_print(b);
    }

  meml_matrix_free (A);
  meml_vector_free (b);

  return (result);
}

#endif DOXYGEN_SHOULD_SKIP_THIS

void qfel_compute_derivatives (const VECTOR * u1, MESH * mesh,
			       MEML_INT node, MEML_FLOAT * u1_xx,
			       MEML_FLOAT * u1_yy, MEML_FLOAT * u1_xy,
			       MEML_FLOAT * u1_x, MEML_FLOAT * u1_y)
{
  INDEXARRAY *liste;
  VECTOR *derivatives;

  /* indizes der nachbarn holen */

  liste = q_hole_indizes (mesh, node);

  derivatives = apply_taylor (node, liste, mesh, u1);

  (*u1_x) = derivatives->data[0];
  (*u1_y) = derivatives->data[1];
  (*u1_xx) = derivatives->data[2];
  (*u1_xy) = derivatives->data[3];
  (*u1_yy) = derivatives->data[4];

  meml_vector_free (derivatives);
  meml_indexarray_free (liste);
}

void qfel_grad_by_taylor (const VECTOR * v, MESH * mesh, VECTOR * v_x,
			  VECTOR * v_y)
{
  MEML_INT i;
  MEML_FLOAT u1_xx, u1_yy, u1_xy, u1_x, u1_y;

  for (i = 0; i < v->dim; i++)
    {
      qfel_compute_derivatives (v, mesh, i, &u1_xx, &u1_yy, &u1_xy, &u1_x,
				&u1_y);
      v_x->data[i] = u1_x;
      v_y->data[i] = u1_y;
    }
}

void qfel_derivatives_by_taylor (const VECTOR * v, MESH * mesh,
				 VECTOR * v_x, VECTOR * v_y,
				 VECTOR * v_xx, VECTOR * v_xy, VECTOR * v_yy)
{
  MEML_INT i;
  MEML_FLOAT u1_xx, u1_yy, u1_xy, u1_x, u1_y;

  for (i = 0; i < v->dim; i++)
    {
      qfel_compute_derivatives (v, mesh, i, &u1_xx, &u1_yy, &u1_xy, &u1_x,
				&u1_y);
      v_x->data[i] = u1_x;
      v_y->data[i] = u1_y;
      v_xx->data[i] = u1_xx;
      v_yy->data[i] = u1_yy;
      v_xy->data[i] = u1_xy;
    }
}

VECTOR * qfel_adapt_lfel_data(VECTOR * ldata, MESH * mesh)
{
    MEML_INT i;
    MEML_FLOAT z1,z2;
    VECTOR * qdata;
    
    qdata = meml_vector_copy(ldata);
    vecl_vector_resize_f (qdata, mesh->number_of_points + mesh->number_of_adof);
    
    for (i = 0; i < mesh->number_of_edges; i++)
    {
	
	z1 = ldata->data[mesh->edges[i][0] - 1];
	z2 = ldata->data[mesh->edges[i][1] - 1];
	
	qdata->data[i+ mesh->number_of_points] = (z1 + z2)/2;
    }

    return (qdata);
}