lfel.c

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

  meml_vector_free(r);

  if ( level-1 == lastLevel)
    {
#ifdef __HAVE_UMFPACK	
	xDown = disl_umfpack (A[lastLevel],rDown);
#else
	xDown = meml_vector_new(rDown->dim);
      if (A[lastLevel]->symmetric=='y')
	itsl_minres (A[lastLevel], rDown,xDown,1e-12);
      else
	itsl_gmres_jacobi (A[lastLevel], rDown,xDown,1e-12,50);
#endif

    }
  else
    {
      xDown = meml_vector_new(rDown->dim);
      lfel_multigrid_function(A, rDown, xDown,mesh,level-1, lastLevel,k1,k2);
    }
  
  meml_vector_free(rDown);

  x = meml_matrix_t_vector_mul (mesh[level]->restriction,xDown);
  
  meml_vector_free(xDown);


  meml_vector_add_f(x,guess);
  
  meml_vector_free(x);

  for (i=0;i<k2;i++)
    glaetter(A[level], b, guess);
  
}

double lfel_multigrid_solver(ME_MATRIX ** A, VECTOR * b, VECTOR * guess, 
			     MESH ** mesh, MEML_INT level, MEML_INT lastLevel,
			     MEML_INT k1,MEML_INT k2, MEML_FLOAT eps, 
			     MEML_INT maxiterations)
{
  double rneu=eps+1;
  double ralt = 0;
  VECTOR * r, *temp;
  int durchlauf = 0,i;

  if ( level <= lastLevel)
    {
#ifdef __HAVE_UMFPACK	
	temp = disl_umfpack (A[level],b);
	meml_vector_copy_f (temp,guess);
	meml_vector_free(temp);
	return(0);
#else
      if (A[lastLevel]->symmetric=='y')
	itsl_minres (A[level], b,guess,eps);
      else
	itsl_gmres_jacobi (A[level], b,guess,1e-12,eps);
#endif
      return(1e-12);
    }


  r = meml_matrix_vector_mul(A[level],b);
  meml_vector_add_ff(-1,b,r);
  rneu = meml_vector_norm(r);
  meml_vector_free(r);
  ralt = rneu+1;

  if (MEML_VERBOSE)
    {
      printf ("ML : r%d = %e \n", 0,rneu);
    }
  
  while ( (rneu>eps) && (rneu <ralt) && (maxiterations>durchlauf) )
    {
      durchlauf++;
      temp = meml_vector_copy(guess);
      lfel_multigrid_function(A, b, temp, mesh, level, lastLevel,k1,k2);      

      r = meml_matrix_vector_mul(A[level],temp);
      meml_vector_add_ff(-1,b,r);
      ralt = rneu;
      rneu = meml_vector_norm(r);
      meml_vector_free(r);
      if (MEML_VERBOSE)
	{
	  printf ("ML : r%d = %e \n", durchlauf,rneu);
	}

      if ( rneu >= ralt )
	{
	  printf("Multigrid divergiert bei %e nach %d Schritten \n",
		 rneu,durchlauf);
	  printf("liefere Naeherung mit %e \n",ralt);
	  meml_vector_free(temp);
	
#ifdef __HAVE_UMFPACK	
	  temp = disl_umfpack (A[level],b);
	  meml_vector_copy_f (temp,guess);
	  meml_vector_free(temp);
	  return(0);
#else
	  if (A[lastLevel]->symmetric=='y')
	    itsl_minres (A[level], b,guess,eps);
	  else
	    itsl_gmres_jacobi (A[level], b,guess,1e-12,eps);
#endif
	}
      else
	{
	  for (i=0;i<guess->dim;i++)
	    guess->data[i] = temp->data[i];
	  meml_vector_free(temp);
	}
    }

  return (rneu);
}




/**
 * \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 lfel_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);
    }
}

/**
 * \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 *lfel_apply_function_to_vector (MESH * mesh, xyt_function f,
				       MEML_FLOAT t)
{
  VECTOR *v;

  v = vecl_vector_new (mesh->number_of_points);

  lfel_apply_function_to_vector_f (v, mesh, f, t);

  return (v);
}

/**
 * \brief  calculates the Neumann boundary condition form a vector field 
 *         and the mesh data
 *
 * \param v vector field at all nodes. 
 *        Only the values at the boundary nodes will be used.
 * \param mesh the mesh 
 * 
 */
VECTOR *lfel_calculat_neumann_boundary_condition (VECTOR ** v, MESH * mesh)
{
  VECTOR *g;

  g = meml_vector_new (mesh->number_of_points);

  lfel_calculat_neumann_boundary_condition_f (v, mesh, g);

  return (g);
}

/**
 * \brief  calculates the Neumann boundary condition form a vector field 
 *         and the mesh data
 *
 * \param v vector field at all nodes. 
 *          Only the values at the boundary nodes will be used.
 * \param mesh the mesh 
 * \param g vector which contains the Neumann boundary condition on exit. 
 *        g has the dimension number_of_points
 * 
 */
void lfel_calculat_neumann_boundary_condition_f (VECTOR ** v, MESH * mesh,
						 VECTOR * g)
{
  MEML_INT i, p1, p2, which_edge;
  MEML_FLOAT temp;

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

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

      temp =
	mesh->normal_vector[i]->data[0] * v[p1 - 1]->data[0] +
	mesh->normal_vector[i]->data[1] * v[p1 - 1]->data[1];
      g->data[p1 - 1] += temp;

      temp =
	mesh->normal_vector[i]->data[0] * v[p2 - 1]->data[0] +
	mesh->normal_vector[i]->data[1] * v[p2 - 1]->data[1];
      g->data[p2 - 1] += temp;
    }

  meml_vector_scaling_f (0.5, g);
}


/**
 * \brief 
 * applies the Neumann boundary compatibility condition to the right side 
 *
 *  g must be a vector of number_of_points size. 
 *  Set for every boundary point i the value of g[i].
 *  
 */
void lfel_set_compatibility_condition_to_right_side (VECTOR * b,
						     const VECTOR * g,
						     MESH * mesh)
{
  MEML_FLOAT S, mu_omega, c;
  MEML_INT number_of_points = b->dim;
  VECTOR *Ones;
  MEML_INT i, j, aktueller_knoten, zaehler, nachbarn[2];
  MEML_FLOAT x[3], y[3], h[2], S_rand = 0;

  if (g!=NULL)
    {
      for (i = 0; i < mesh->number_of_boundary_points; i++)
	{
	  aktueller_knoten = mesh->boundary->data[i];
	  zaehler = 0;

	  /* finde die rand-nachbarn des aktuellen knotens */
	  for (j = 0; j < mesh->number_of_neighbours[aktueller_knoten - 1]; j++)
	    if (meml_element_of_indexarray
		(mesh->boundary, mesh->neighbours[aktueller_knoten - 1]->data[j]))
	      {
		nachbarn[zaehler] =
		  mesh->neighbours[aktueller_knoten - 1]->data[j];
		zaehler++;
	      }
	  /* bestimme den Abstand */
	  for (j = 0; j < 2; j++)
	    {
	      x[j] = mesh->points->data[2 * (nachbarn[j] - 1)];
	      y[j] = mesh->points->data[2 * (nachbarn[j] - 1) + 1];
	    }
	  
	  x[2] = mesh->points->data[2 * (aktueller_knoten - 1)];
	  y[2] = mesh->points->data[2 * (aktueller_knoten - 1) + 1];
	  
	  for (j = 0; j < 2; j++)
	    {
	      h[j] = sqrt (pow (x[j] - x[2], 2) + pow (y[j] - y[2], 2));
	      /* Simpson-Regel */
	      S_rand = S_rand + h[j] * (0.5 * g->data[aktueller_knoten - 1] +
				    0.5 * g->data[nachbarn[j] - 1]);
	    }
	}
      
      /* da jedes stueck doppel gezaehlt wurde */
      S_rand = S_rand / 2;
    }
  else
    {
      S_rand = 0;
    }

  /* compute the integral of f over omega */
  S = lfel_integrating (b, mesh);

  /* compute \mu(omega) */
  Ones = meml_vector_new_ones (b->dim);
  mu_omega = lfel_integrating (Ones, mesh);
  meml_vector_free (Ones);

  c = (S + S_rand) / mu_omega;

  for (i = 0; i < number_of_points; i++)
    b->data[i] = b->data[i] - c;
}


/**
 * \brief applies the neumann boundary conditions to the right side 
 *
 *  g must be a vector of number_of_points size. 
 *  Set for every boundary point i the value of g[i].
 *  
 */
void lfel_set_neumann_cond_to_right_side (VECTOR * b, const VECTOR * g,
					  MESH * mesh)
{
  MEML_INT i, j, aktueller_knoten, zaehler, nachbarn[2];
  MEML_FLOAT x[3], y[3], h[2], integral_wert[2];

  for (i = 0; i < mesh->number_of_boundary_points; i++)
    {
      aktueller_knoten = mesh->boundary->data[i];
      zaehler = 0;

      /* finde die rand-nachbarn des aktuellen knotens */
      for (j = 0; j < mesh->number_of_neighbours[aktueller_knoten - 1]; j++)
	{
	  if (meml_element_of_indexarray (mesh->boundary,
					  mesh->
					  neighbours[aktueller_knoten -
						     1]->data[j]))
	    {
	      nachbarn[zaehler] =
		mesh->neighbours[aktueller_knoten - 1]->data[j];
	      zaehler++;
	    }
	}
      /* bestimme den Abstand */
      for (j = 0; j < 2; j++)
	{
	  x[j] = mesh->points->data[2 * (nachbarn[j] - 1)];
	  y[j] = mesh->points->data[2 * (nachbarn[j] - 1) + 1];
	}

      x[2] = mesh->points->data[2 * (aktueller_knoten - 1)];
      y[2] = mesh->points->data[2 * (aktueller_knoten - 1) + 1];

      for (j = 0; j < 2; j++)
	{
	  h[j] = sqrt (pow (x[j] - x[2], 2) + pow (y[j] - y[2], 2));
	  /* Simpson-Regel */
	  integral_wert[j] =
	    h[j] * ((1.0 / 6.0) * g->data[aktueller_knoten - 1] +
		    (2.0 / 3.0) * 0.25 * (g->data[aktueller_knoten - 1] +
					  g->data[nachbarn[j] - 1]));
	  b->data[aktueller_knoten - 1] =
	    b->data[aktueller_knoten - 1] + integral_wert[j];
	}
    }
}

/**
 *  \brief 
 *  Computes the interpolation operator from the coarse 
 *  points P0 to the points P1 
 *
 * The function computes a sparse matrix representation of the interpolation
 * operator from a coarse mesh to a finer mesh.
 *
 * \todo In the monent mesh_high must be the globale 
 * red refinement of mesh_low. A more general interpolation must be build.
 * 
 * \param mesh_low the coarse mesh
 * \param mesh_high the finer mesh
 *
 */
ME_MATRIX *lfel_interpolation (MESH * mesh_low, MESH * mesh_high)
{
  VECTOR *P0 = mesh_low->points;
  VECTOR *P1 = mesh_high->points;
  INDEXARRAY *triangles = mesh_high->triangles;

  MEML_INT n0 = P0->dim / 2;
  MEML_INT n1 = P1->dim / 2;
  MEML_INT nt = 0.75 * triangles->dim / 3;
  TRIANGLE the_triangle;
  MEML_INT i;
  ME_MATRIX *I;
  ME_MATRIX *temp_matrix;
  INDEXARRAY *T1;
  INDEXARRAY *T2;
  INDEXARRAY *T3;
  VECTOR *Ones;

  T1 = meml_indexarray_new (nt);
  T2 = meml_indexarray_new (nt);
  T3 = meml_indexarray_new (nt);

  /* getting the neighborhood */
  for (i = 0; i < nt; i++)
    {
      the_triangle = smfl_get_triangle (i, triangles->data);
      meml_indexarray_element_set (T1, i, the_triangle.p[0] - 1);
      meml_indexarray_element_set (T2, i, the_triangle.p[1] - 1);
      meml_indexarray_element_set (T3, i, the_triangle.p[2] - 1);
    }

  Ones = meml_vector_new_ones (nt);
  meml_vector_scaling_f (0.5, Ones);

  I = meml_sparse (T2, T1, Ones);
  temp_matrix = meml_sparse (T3, T1, Ones);

  meml_matrix_resize (&I, n1, n0);
  meml_matrix_resize (&temp_matrix, n1, n0);

  meml_matrix_add_f (temp_matrix, I);
  meml_matrix_ceil (I);

  meml_matrix_scaling_f (0.5, I);

  meml_matrix_free (temp_matrix);

  temp_matrix = meml_matrix_new_eye (LS, n1, n0);

  meml_matrix_add_f (temp_matrix, I);

  meml_vector_free (Ones);
  meml_indexarray_free (T1);
  meml_indexarray_free (T2);
  meml_indexarray_free (T3);
  meml_matrix_free (temp_matrix);

  return (I);
}


/**
 * \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 
 * \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 lfel_grad (const VECTOR * v, MESH * mesh, VECTOR * v_x, VECTOR * v_y)
{
  TRANSFORMATION *transformation = mesh->transformation;
  INDEXARRAY *triangles = mesh->triangles;

  VECTOR **gradient;
  VECTOR *c;
  VECTOR *b1;
  VECTOR *b2;
  MEML_INT number_of_triangles = (triangles->dim) / 3;
  MEML_INT number_of_points = v->dim;
  MEML_INT i, j;
  MEML_INT *zaehler;
  MEML_FLOAT gradx, grady;
  TRIANGLE the_triangle;
  MEML_INT *triangle_node = triangles->data;

  gradient = (VECTOR **) calloc (number_of_points, sizeof (VECTOR *));
  zaehler = (MEML_INT *) calloc (number_of_points, sizeof (MEML_INT));

  for (i = 0; i < number_of_points; i++)
    {
      gradient[i] = meml_vector_new (2);
    }

  c = meml_vector_new (3);
  b1 = meml_vector_new (3);
  b2 = meml_vector_new (3);

  for (i = 0; i < number_of_triangles; i++)
    {
      /* b1xb2=c */
      /* konstruiere b1 und b2 */
      meml_vector_element_set_f (b1, 0, transformation[i].B[0][0]);
      meml_vector_element_set_f (b1, 1, transformation[i].B[1][0]);

      meml_vector_element_set_f (b2, 0, transformation[i].B[0][1]);
      meml_vector_element_set_f (b2, 1, transformation[i].B[1][1]);

      /* dritte komponente berechnen */
      the_triangle = smfl_get_triangle (i, triangle_node);

      b1->data[2] = v->data[the_triangle.p[1] - 1] -
	v->data[the_triangle.p[0] - 1];
      b2->data[2] = v->data[the_triangle.p[2] - 1] -
	v->data[the_triangle.p[0] - 1];

      /* kreuprodukt berechnen */
      lfel_kreuzprodukt (b1, b2, c);

      if (c->data[2] != 0)
	{
	  gradx = -(c->data[0] / c->data[2]);
	  grady = -(c->data[1] / c->data[2]);

	  for (j = 0; j < 3; j++)
	    {
	      (gradient[the_triangle.p[j] - 1])->data[0] += gradx;