lfel.c

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

	      (gradient[the_triangle.p[j] - 1])->data[1] += grady;
	      zaehler[the_triangle.p[j] - 1]++;
	    }
	}
      else
	{
	  fprintf (stderr, "Error in lfel_grad: div with zero!\n");
	  exit (EXIT_FAILURE);
	}

    }

  meml_vector_free (c);
  meml_vector_free (b1);
  meml_vector_free (b2);

  for (i = 0; i < number_of_points; i++)
    {
	(gradient[i])->data[0] = ((gradient[i])->data[0]) / zaehler[i];
	(gradient[i])->data[1] = ((gradient[i])->data[1]) / zaehler[i];
    }

  free (zaehler);

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

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

  free (gradient);
}

void lfel_add_steamline_tearms(MEML_FLOAT eps, VECTOR *c1 , VECTOR *c2, 
			       MEML_FLOAT r, 
			       MESH * mesh, ME_MATRIX * A)
{
  MEML_FLOAT kmax=0, k[2], kx=0, ky=0;
  MEML_INT i,x, number_of_triangles = mesh->number_of_triangles;
  TRIANGLE the_triangle;
  MEML_FLOAT delta,Pe;
  VECTOR *cc1, *cc2, *cr;
  VECTOR *c1x, *c1y, *c2x, *c2y;

  if (c1->dim != c2->dim)
    {
      fprintf(stderr,"lfel_add_steamline_tearms : 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_add_steamline_tearms : size of c1 can be 1, "
	      "the number of points or the number of triangles.\n");
      return;
    }

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

  /* es gibt keine konvection */
  if (kmax==0)
    {
      return;
    }
  
  cc1 = meml_vector_new(number_of_triangles);
  cc2 = meml_vector_new(number_of_triangles);
  cr  = meml_vector_new(number_of_triangles);

  c1x = meml_vector_new(mesh->number_of_points);
  c2x = meml_vector_new(mesh->number_of_points);
  c1y = meml_vector_new(mesh->number_of_points);
  c2y = meml_vector_new(mesh->number_of_points);

  if ( c1->dim == mesh->number_of_points)
    {
      lfel_grad_by_taylor (c1, mesh,c1x,c1y);
      lfel_grad_by_taylor (c2, mesh,c2x,c2y);
    }

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

      if (c1->dim == mesh->number_of_points) 
	{
	  k[0]=0;
	  k[1]=0;
	  kx = 0;
	  ky = 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];
	      kx = c1x->data[the_triangle.p[x] - 1];
	      ky = c2y->data[the_triangle.p[x] - 1];
	    }
	  
	  k[0] = k[0]/3;
	  k[1] = k[1]/3;

	  kx = kx/3;
	  ky = ky/3;

	}
      else if (c1->dim == mesh->number_of_triangles) 
	{
	  fprintf(stderr,"lfel_add_steamline_tearms : dim of c1 = "
		  "number of triangles is not complete until now\n");

	  k[0] = c1->data[i];
	  k[1] = c2->data[i];
	}
      else /* (c1->dim != 1) ) */ 
	{
	  k[0] = c1->data[0];
	  k[1] = c2->data[0];
	  
	  kx = 0;
	  ky = 0;
	}
      
      Pe = (meml_max(fabs(k[0]),fabs(k[1]))) * mesh->ht[i] / ( 2*eps);
      
      if (Pe>1)
	{
	  delta = mesh->ht[i] * f_delta1(r,kmax,mesh->hmax,eps);
	}
      else
	{
	  delta = 0;
	}
      
      cc1->data[i] = - r * delta * k[0];
      cc2->data[i] = - r * delta * k[1];

      cr->data[i] = -delta * (kx + ky) * r;
    }

  lfel_assemb_convection(cc1,cc2,A,mesh);
  
  if ( c1->dim != 1)
    if ( r != 0)
      lfel_assemb_reaction_f(cr,A,mesh);

  meml_vector_free(cc1);
  meml_vector_free(cc2);
  meml_vector_free(cr);

  meml_vector_free(c1x);
  meml_vector_free(c1y);
  meml_vector_free(c2y);
  meml_vector_free(c2x);
}

void lfel_assemb_pde(MEML_FLOAT eps, VECTOR *c1 , VECTOR *c2, MEML_FLOAT r, 
		     xyt_function f1, MEML_FLOAT t, VECTOR * f2, 
		     MESH * mesh, ME_MATRIX ** A, VECTOR ** b)
{
  ME_MATRIX * B;
  VECTOR *b1, *b2, *bs1, *bs2;

  B=meml_matrix_new(LS,mesh->number_of_points,mesh->number_of_points);

  if (eps != 0)
    {
      lfel_assemb_diffusion(eps,B,mesh);
    }

  if (r!=0)
    {
      lfel_assemb_reaction(r,B,mesh);
    }

  if ( c1 != NULL)
    {
      lfel_assemb_convection(c1,c2,B,mesh);
      lfel_assemb_steamline_diffusion(eps,c1,c2,r,B,mesh);
      lfel_add_reaction_steamline_term (B,eps,c1 ,c2,r,mesh);
    }

  if (A!=NULL)
    if ( (eps != 0) || (r!=0) || ( c1 != NULL) )
      (*A)=meml_matrix_convert(CS,B);
  
  meml_matrix_free(B);

  if (f1 != NULL)
    {
      b1=lfel_assemb_load_vector(f1,t,mesh);
      if ( c1 !=NULL)
	{
	  bs1 = lfel_assemb_steamline_diffusion_rightside 
	      (eps,c1,c2,r,f1,t,mesh);
	}
      else
	bs1=meml_vector_new(mesh->number_of_points);
    }
  else
    {
      b1=meml_vector_new(mesh->number_of_points);
      bs1=meml_vector_new(mesh->number_of_points);
    }

  if ( f2 != NULL)
    {
      b2= lfel_assemb_right_side (f2, mesh);
      if ( c1 !=NULL)
	{
	  bs2 = lfel_assemb_steamline_diffusion_rightside_v 
	      (eps,c1,c2,r,f2,mesh);
	}
      else
	bs2=meml_vector_new(mesh->number_of_points);
    }
  else
    {
      b2=meml_vector_new(mesh->number_of_points);
      bs2=meml_vector_new(mesh->number_of_points);
    }
  
  if (b!=NULL)
    {
      (*b)=meml_vector_add(b1,b2);
      meml_vector_add_ff(1,bs1,(*b));
      meml_vector_add_ff(1,bs2,(*b));
    }

  meml_vector_free(b1);
  meml_vector_free(b2); 
  meml_vector_free(bs1);
  meml_vector_free(bs2);
}

void lfel_assemb_steamline_diffusion(MEML_FLOAT eps, VECTOR * c1, VECTOR * c2,
				     MEML_FLOAT r, ME_MATRIX * A, 
				     const MESH * mesh)
{
    lfel_assemb_diffusion_ff(eps, c1, c2, r, A,mesh,'y');
}

void lfel_assemb_diffusion_ff(MEML_FLOAT eps, VECTOR * c1, VECTOR * c2,
			      MEML_FLOAT r, ME_MATRIX * A, 
			      const MESH * mesh, const char UseAsSteamline)
{
  TRANSFORMATION *transformation = mesh->transformation;
  MEML_FLOAT b[3][2],k[2],E[3][3],BT[2][2],det,kmax=0;
  MEML_INT i,x,y, number_of_triangles = mesh->number_of_triangles;
  TRIANGLE the_triangle;
  MEML_FLOAT delta,Pe;

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

  /* es gibt keine konvection */
  if (kmax==0)
    {
      return;
    }
  
  for (i=0;i<number_of_triangles;i++)
    {
      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 ( UseAsSteamline == 'n')
      {
	  Pe = 2;
	  delta = 1;
      }
      else
      {
	  Pe = meml_max(k[0],k[1]) * mesh->ht[i] / ( 2*eps);
	  delta = mesh->ht[i] * f_delta1(r,kmax,mesh->hmax,eps);
      }
      
      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]' */

	  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]) *
				(k[0]*b[y][0]+k[1]*b[y][1]) )
			/(2*fabs(transformation[i].det));
		}
	    }
	  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, delta*E[x][y]);
		}
	    }
	}
    }

}


void lfel_add_reaction_steamline_term (ME_MATRIX * A,
				       MEML_FLOAT eps, VECTOR *c1 , VECTOR *c2 ,
				       MEML_FLOAT r, const MESH * mesh)
{
  TRANSFORMATION *transformation = mesh->transformation;
  MEML_FLOAT b[3][2],k[2],BT[2][2],det;
  MEML_INT i,x,y, number_of_triangles = mesh->number_of_triangles;
  TRIANGLE the_triangle;
  MEML_FLOAT delta,kmax=0,Pe;
  
  if (c1->dim != c2->dim)
    {
      fprintf(stderr,"lfel_add_reaction_steamline_term : "
	      "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_add_reaction_steamline_term : size of c1 can be 1, "
	      "the number of points or the number of triangles.\n");
      exit(EXIT_FAILURE);
    }
  
  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;
    }

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

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

      Pe = meml_max(k[0],k[1]) * mesh->ht[i] / ( 2*eps);

      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 < 3; x++)
	    for (y = 0; y < 3; y++)
	      {
		meml_matrix_element_add (A, the_triangle.p[x] - 1,
					 the_triangle.p[y] - 1, 
					 r* delta * (k[0]*b[x][0]+k[1]*b[x][1])/6.0);
	      }
	}
    }
}

VECTOR * lfel_assemb_steamline_diffusion_rightside_v(MEML_FLOAT eps,VECTOR *c1, 
						     VECTOR *c2 , MEML_FLOAT r,
						     VECTOR *f, const MESH *mesh)
{
  return ( 
	  lfel_assemb_steamline_diffusion_rightside_vf(eps,c1, c2 , r, f,mesh,'n')