lfel.c

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

	  );
}


VECTOR * lfel_assemb_steamline_diffusion_rightside_vf(MEML_FLOAT eps,VECTOR *c1, 
						      VECTOR *c2 , MEML_FLOAT r,
						      VECTOR *f,const MESH *mesh,
						      const char IgnorePe)
{
  TRANSFORMATION *transformation = mesh->transformation;
  MEML_FLOAT b[3][2],k[2],BT[2][2],det, f_temp;
  MEML_INT i,x, number_of_triangles = mesh->number_of_triangles;
  TRIANGLE the_triangle;
  VECTOR *temp_vector;
  MEML_FLOAT delta,kmax=0,Pe;

  if (c1->dim != c2->dim)
    {
      fprintf(stderr,"lfel_assemb_convection : "
	      "size of c1 and c2 muss agree!\n");
      exit(EXIT_FAILURE);
    }

  if ( (c1->dim != mesh->number_of_points) && 
       (c1->dim != mesh->number_of_triangles) && 
       (c1->dim != 1) )
    {
      fprintf(stderr,"lfel_assemb_convection : size of c1 can be 1, "
	      "the number of points or the number of triangles.\n");
      exit(EXIT_FAILURE);
    }

  temp_vector = meml_vector_new(mesh->number_of_points);

  for (i=0;i<c1->dim;i++)
    {
      k[0] = fabs(c1->data[i]);
      k[1] = fabs(c2->data[i]);
      
      if (k[0]>kmax)
	kmax = k[0];
      if (k[1]>kmax)
	kmax = k[1];
    }

  /* nullvektor zurueck geben */
  if (kmax==0)
    {
      return(temp_vector);
    }

  for (i=0;i<number_of_triangles;i++)
    {
      if (IgnorePe == 'n')
	delta = mesh->ht[i] * f_delta1(r,kmax,mesh->hmax,eps);
      else
	delta = 1;

      the_triangle = smfl_get_triangle(i, mesh->triangles->data);

      if (c1->dim == mesh->number_of_points) 
	{
	  k[0]=0;
	  k[1]=0;
	  
	  for (x=0;x<3;x++)
	    {
	      k[0] += c1->data[the_triangle.p[x] - 1];
	      k[1] += c2->data[the_triangle.p[x] - 1];
	    }
      
	  k[0] = k[0]/3;
	  k[1] = k[1]/3;
	}
      else if (c1->dim == mesh->number_of_triangles) 
	{
	  k[0] = c1->data[i];
	  k[1] = c2->data[i];
	}
      else /* (c1->dim != 1) ) */ 
	{
	  k[0] = c1->data[0];
	  k[1] = c2->data[0];
	}

      if (IgnorePe == 'n')
	Pe = meml_max(k[0],k[1]) * mesh->ht[i] / ( 2*eps);
      else
	Pe = 42;

      if (Pe>1)
	{
	  det = transformation[i].det/fabs(transformation[i].det);

	  /*
	    B =  
	    ( B00 B01 )
	    ( B10 B11 )
	    
	    B^-T = ( B11 -B01 ) 
	           (-B10  B00 )    * 1/det  
	  */
	  
	  BT[0][0] =  transformation[i].B[1][1]*det;
	  BT[0][1] = -transformation[i].B[1][0]*det;
	  BT[1][0] = -transformation[i].B[0][1]*det;
	  BT[1][1] =  transformation[i].B[0][0]*det;
	  
	  
	  b[0][0] =  -BT[0][0]-BT[0][1];
	  b[0][1] =  -BT[1][0]-BT[1][1];
      
	  b[1][0] =   BT[0][0];
	  b[1][1] =   BT[1][0];
	  
	  b[2][0] =   BT[0][1];
	  b[2][1] =   BT[1][1];
	  
	  /* K = [ k0, k1] * [k0,k1]' */
	  
	  the_triangle = smfl_get_triangle (i, mesh->triangles->data);
	  
	  f_temp= 0;
	  
	  for (x = 0; x < 3; x++)
	    {
	      f_temp += f->data[the_triangle.p[x] - 1];
	    }

	  f_temp=(f_temp/6);
	  
	  for (x = 0; x < 3; x++)
	    {
	      temp_vector->data[the_triangle.p[x] - 1] += delta * 
		f_temp * (k[0]*b[x][0]+k[1]*b[x][1]);
	    }
	}
    }
 
  return (temp_vector);
}

VECTOR * lfel_assemb_right_side_least_squares (VECTOR * c1,
					       VECTOR * c2 , VECTOR * r,
					       xyt_function f1, VECTOR *f2,  
					       MEML_FLOAT time, 
					       const MESH * mesh)
{
    VECTOR *b1 , *b2;

    if ( (r->dim != mesh->number_of_points) && (r->dim != mesh->number_of_triangles) && (r->dim != 1) )
    {
	fprintf(stderr,"lfel_assemb_right_side_least_squares : size of r can be 1, "
		"the number of points or the number of triangles.\n");
	exit(EXIT_FAILURE);
    }

    if (f1 != NULL)
      {
	b1 = lfel_assemb_steamline_diffusion_rightside_f (1, c1, c2 , 1, 
							  f1, time, mesh, 'y');

	
	b2 = lfel_assemb_load_vector_f (f1,time,r, mesh);
	meml_vector_add_f (b2, b1);
	meml_vector_free(b2);
      }
    else
      {
	b1 = meml_vector_new(mesh->number_of_points);
      }


    if (f2 != NULL)
      {
	b2 = lfel_assemb_right_side_f (f2, r, mesh);
	meml_vector_add_f (b2, b1);
	meml_vector_free(b2);
	
	b2 = lfel_assemb_steamline_diffusion_rightside_vf(1,c1,c2,1,f2,mesh, 'y');
	meml_vector_add_f (b2, b1);
	meml_vector_free(b2);
      }

    return(b1);
}

VECTOR * lfel_assemb_steamline_diffusion_rightside (MEML_FLOAT eps, VECTOR * c1,
						    VECTOR * c2 , MEML_FLOAT r,
						    xyt_function f, 
						    MEML_FLOAT time, 
						    const MESH * mesh)
{

    return( 
	lfel_assemb_steamline_diffusion_rightside_f (eps,c1,c2 ,r,f,time,mesh,'n')
	);
    
}


VECTOR * lfel_assemb_steamline_diffusion_rightside_f (MEML_FLOAT eps, VECTOR * c1,
						      VECTOR * c2 , MEML_FLOAT r,
						      xyt_function f, 
						      MEML_FLOAT time, 
						      const MESH * mesh, 
						      const char IgnorePe)
{
  TRANSFORMATION *transformation = mesh->transformation;
  MEML_FLOAT b[3][2],k[2],BT[2][2],det, d[2], f_temp;
  MEML_INT i,x, number_of_triangles = mesh->number_of_triangles;
  TRIANGLE the_triangle;
  VECTOR *points = mesh->points;
  VECTOR *temp_vector;
  MEML_FLOAT f_value[7];
  MEML_FLOAT delta,kmax=0,Pe;

  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
  };

  if (c1->dim != c2->dim)
    {
      fprintf(stderr,"lfel_assemb_convection : "
	      "size of c1 and c2 muss agree!\n");
      exit(EXIT_FAILURE);
    }

  if ( (c1->dim != mesh->number_of_points) && 
       (c1->dim != mesh->number_of_triangles) && 
       (c1->dim != 1) )
    {
      fprintf(stderr,"lfel_assemb_convection : size of c1 can be 1, "
	      "the number of points or the number of triangles.\n");
      exit(EXIT_FAILURE);
    }

  temp_vector = meml_vector_new(mesh->number_of_points);

  for (i=0;i<c1->dim;i++)
    {
      k[0] = fabs(c1->data[i]);
      k[1] = fabs(c2->data[i]);
      
      if (k[0]>kmax)
	kmax = k[0];
      if (k[1]>kmax)
	kmax = k[1];
    }

  /* nullvektor zurueck geben */
  if (kmax==0)
    {
      return(temp_vector);
    }

  for (i=0;i<number_of_triangles;i++)
    {
      if (IgnorePe == 'n')
	  delta = mesh->ht[i] * f_delta1(r,kmax,mesh->hmax,eps);
      else
	  delta = 1;

      the_triangle = smfl_get_triangle(i, mesh->triangles->data);

      if (c1->dim == mesh->number_of_points) 
	{
	  k[0]=0;
	  k[1]=0;
	  
	  for (x=0;x<3;x++)
	    {
	      k[0] += c1->data[the_triangle.p[x] - 1];
	      k[1] += c2->data[the_triangle.p[x] - 1];
	    }
      
	  k[0] = k[0]/3;
	  k[1] = k[1]/3;
	}
      else if (c1->dim == mesh->number_of_triangles) 
	{
	  k[0] = c1->data[i];
	  k[1] = c2->data[i];
	}
      else /* (c1->dim != 1) ) */ 
	{
	  k[0] = c1->data[0];
	  k[1] = c2->data[0];
	}

      if (IgnorePe == 'n')
	  Pe = meml_max(k[0],k[1]) * mesh->ht[i] / ( 2*eps);
      else
	  Pe = 42;

      if (Pe>1)
	{
	  det = transformation[i].det/fabs(transformation[i].det);

	  /*
	    B =  
	    ( B00 B01 )
	    ( B10 B11 )
	    
	    B^-T = ( B11 -B01 ) 
	    (-B10  B00 )    * 1/det  
	  */
	  
	  BT[0][0] =  transformation[i].B[1][1]*det;
	  BT[0][1] = -transformation[i].B[1][0]*det;
	  BT[1][0] = -transformation[i].B[0][1]*det;
	  BT[1][1] =  transformation[i].B[0][0]*det;
	  
	  
	  b[0][0] =  -BT[0][0]-BT[0][1];
	  b[0][1] =  -BT[1][0]-BT[1][1];
      
	  b[1][0] =   BT[0][0];
	  b[1][1] =   BT[1][0];
	  
	  b[2][0] =   BT[0][1];
	  b[2][1] =   BT[1][1];
	  
	  /* K = [ k0, k1] * [k0,k1]' */
	  
	  the_triangle = smfl_get_triangle (i, mesh->triangles->data);
	  
	  for (x = 0; x < 7; x++)
	    {
	      d[0] = points->data[(the_triangle.p[0] - 1) * 2];
	      d[1] = points->data[(the_triangle.p[0] - 1) * 2 + 1];
	      
	      f_value[x] = (*f) (transformation[i].B[0][0] * gaus_x[x] +
				 transformation[i].B[0][1] * gaus_y[x] + d[0],
				 transformation[i].B[1][0] * gaus_x[x] +
				 transformation[i].B[1][1] * gaus_y[x] + d[1],
				 time);
	    }
	  
	  f_temp= 0;
	  
	  for (x = 0; x < 7; x++)
	    {
	      f_temp += gaus_w[x] * f_value[x];
	    }
	  
	  for (x = 0; x < 3; x++)
	    {
	      temp_vector->data[the_triangle.p[x] - 1] += delta * 
		f_temp * (k[0]*b[x][0]+k[1]*b[x][1]);
	    }
	}
    }
 
  return (temp_vector);
}

void lfel_assemb_convection(VECTOR *c1 , VECTOR *c2, 
			    ME_MATRIX * A, const MESH * mesh)
{
  TRANSFORMATION *transformation = mesh->transformation;
  MEML_INT i,x,y, number_of_triangles = mesh->number_of_triangles;
  MEML_FLOAT k[2], b[3][2], det,BT[2][2],E[3][3];
  TRIANGLE the_triangle;

  if (c1->dim != c2->dim)
    {
      fprintf(stderr,"lfel_assemb_convection : size of c1 and c2 muss agree!\n");
      return;
    }

  if ( (c1->dim != mesh->number_of_points) && 
       (c1->dim != mesh->number_of_triangles) && 
       (c1->dim != 1) )
    {
      fprintf(stderr,"lfel_assemb_convection : size of c1 can be 1, "
	      "the number of points or the number of triangles.\n");
      return;
    }

  for (i=0;i<number_of_triangles;i++)
    {
      the_triangle = smfl_get_triangle(i, mesh->triangles->data);
      
      det = transformation[i].det/fabs(transformation[i].det);
      
      /*
	B =  
               ( B00 B01 )
	       ( B10 B11 )

	B^-T = ( B11 -B01 ) 
               (-B10  B00 )    * 1/det  
      */

      BT[0][0] =  transformation[i].B[1][1]*det;
      BT[0][1] = -transformation[i].B[1][0]*det;
      BT[1][0] = -transformation[i].B[0][1]*det;
      BT[1][1] =  transformation[i].B[0][0]*det;


      b[0][0] =  -BT[0][0]-BT[0][1];
      b[0][1] =  -BT[1][0]-BT[1][1];

      b[1][0] =   BT[0][0];
      b[1][1] =   BT[1][0];

      b[2][0] =   BT[0][1];
      b[2][1] =   BT[1][1];

      /*
      printf("B=[%e %e ; %e %e]\n",transformation[i].B[0][0],transformation[i].B[0][1]
	     ,transformation[i].B[1][0],transformation[i].B[1][1]);

      printf("BT=[%e %e ; %e %e]\n",BT[0][0],BT[0][1],BT[1][0],BT[1][1]);
      */

      if (c1->dim == mesh->number_of_points) 
	{
	  k[0]=0;
	  k[1]=0;
	  
	  for (x=0;x<3;x++)
	    {
	      k[0] += c1->data[the_triangle.p[x] - 1];
	      k[1] += c2->data[the_triangle.p[x] - 1];
	    }
      
	  k[0] = k[0]/3;
	  k[1] = k[1]/3;
	}
      else if (c1->dim == mesh->number_of_triangles) 
	{
	  k[0] = c1->data[i];
	  k[1] = c2->data[i];
	}
      else /* (c1->dim != 1) ) */ 
	{
	  k[0] = c1->data[0];
	  k[1] = c2->data[0];
	}
      /*
      printf("k auf dreieck no %d = (%e , %e) \n",i,k[0],k[1]);
      */

      for(y=0;y<3;y++)
	{
	  for (x = 0; x < 3; x++)
	    {
	      E[y][x] = (k[0]*b[x][0]+k[1]*b[x][1])/6;
	    }
	}

      for (x = 0; x < 3; x++)
	{
	  for (y = 0; y < 3; y++)
	    {
	      /* wegen C <-> Matlab -1 */
	      meml_matrix_element_add_f(A, the_triangle.p[x] - 1,
					the_triangle.p[y] - 1, E[x][y]);
	    }
	}
    }
}

/**
 *   \brief 
 *   \~german 
 *   Assembliert die Steifigkeitsmatrix/den Diffusionsanteil
 *   fuer das Galerkinverfahren
 *   \~english 
 *   Does the assembly of the global stiffness matrix 
 *
 *   \~german 
 *   Die Matrix A muss vor dem Aufruf dieser Funktion erzeugt worden sein.
 *   Bei groesseren Gitter sollte A eine Sparse-Matrix-Typ sein.
 *   \~english
 *   The matrix A got be initialized before calling lfel_assemb_diffusion.
 *   On a finer mesh you should use a sparce type for A.
 *  
 *   \~german 
 *  \param eps Diffusionsparameter der partiellen Differentialgleichung
 *  \param A die Matrix zu der die Steifigkeitsmatrix addiert werden soll
 */
int lfel_assemb_diffusion (const MEML_FLOAT eps, ME_MATRIX * A, MESH * mesh)
{
  VECTOR * epsarray;
  int r;

  epsarray = meml_vector_new(1);

  epsarray->data[0]=eps;

  r=lfel_assemb_diffusion_f (epsarray, A, mesh);

  meml_vector_free(epsarray);
  
  return(r);
}

/**
 *   \brief 
 *   \~german 
 *   Assembliert die Steifigkeitsmatrix/den Diffusionsanteil
 *   fuer das Galerkinverfahren
 *   \~english 
 *   Does the assembly of the global stiffness matrix 
 *
 *   \~german 
 *   Die Matrix A muss vor dem Aufruf dieser Funktion erzeugt worden sein.
 *   Bei groesseren Gitter sollte A eine Sparse-Matrix-Typ sein.
 *   \~english
 *   The matrix A got be initialized before calling lfel_assemb_diffusion.
 *   On a finer mesh you should use a sparce type for A.
 *  
 *   \~german 
 *  \param eps Diffusionsparameter der partiellen Differentialgleichung
 *  \param A die Matrix zu der die Steifigkeitsmatrix addiert werden soll