mesh.c

OpenFVM-v1.1 open source cfd code

				 nodes[faces[nbfaces].node[3]]);

	      // Calculate centroid   
	      faces[nbfaces].cface =
		GeoCalcCentroid4 (nodes[faces[nbfaces].node[0]],
				  nodes[faces[nbfaces].node[1]],
				  nodes[faces[nbfaces].node[2]],
				  nodes[faces[nbfaces].node[3]]);

	      // Calculate normal of the face                         
	      faces[nbfaces].n = GeoCalcNormal (nodes[faces[nbfaces].node[0]],
						nodes[faces[nbfaces].node[1]],
						nodes[faces[nbfaces].
						      node[2]]);

	      d =
		GeoSubVectorVector (faces[nbfaces].cface,
				    elements[element].celement);

	      // Normal should point to neighbour
	      // Flip normal if necessary
	      if (GeoDotVectorVector (faces[nbfaces].n, d) < 0.0)
		{
		  faces[nbfaces].n.x *= -1;
		  faces[nbfaces].n.y *= -1;
		  faces[nbfaces].n.z *= -1;
		}

	      faces[nbfaces].element = element;
	      faces[nbfaces].physreg = -1;
	      faces[nbfaces].elemreg = -1;
              faces[nbfaces].partition = -1;
	      faces[nbfaces].pair = -1;
	      faces[nbfaces].bc = NONE;

	      nbfaces++;

	    }
	}

      if (elements[element].type == PRISM)
	{

	  elements[element].nbfaces = 5;

	  // Allocate memory
	  elements[element].face =
	    calloc (elements[element].nbfaces, sizeof (int));

	  for (j = 0; j < 2; j++)
	    {

	      elements[element].face[j] = nbfaces;

	      faces[nbfaces].type = TRIANGLE;
	      faces[nbfaces].nbnodes = 3;

	      switch (j)
		{
		case 0:

		  // Allocate memory
		  faces[nbfaces].node =
		    calloc (faces[nbfaces].nbnodes, sizeof (int));

		  faces[nbfaces].node[0] = elements[element].node[0];
		  faces[nbfaces].node[1] = elements[element].node[1];
		  faces[nbfaces].node[2] = elements[element].node[2];

		  break;

		case 1:

		  // Allocate memory
		  faces[nbfaces].node =
		    calloc (faces[nbfaces].nbnodes, sizeof (int));

		  faces[nbfaces].node[0] = elements[element].node[3];
		  faces[nbfaces].node[1] = elements[element].node[5];
		  faces[nbfaces].node[2] = elements[element].node[4];

		  break;

		}

	      // Calculate area
	      faces[nbfaces].Aj =
		GeoCalcTriArea (nodes[faces[nbfaces].node[0]],
				nodes[faces[nbfaces].node[1]],
				nodes[faces[nbfaces].node[2]]);

	      // Calculate centroid   
	      faces[nbfaces].cface =
		GeoCalcCentroid3 (nodes[faces[nbfaces].node[0]],
				  nodes[faces[nbfaces].node[1]],
				  nodes[faces[nbfaces].node[2]]);

	      // Calculate normal of the face                         
	      faces[nbfaces].n = GeoCalcNormal (nodes[faces[nbfaces].node[0]],
						nodes[faces[nbfaces].node[1]],
						nodes[faces[nbfaces].
						      node[2]]);

	      d =
		GeoSubVectorVector (faces[nbfaces].cface,
				    elements[element].celement);

	      // Normal should point to neighbour
	      // Flip normal if necessary
	      if (GeoDotVectorVector (faces[nbfaces].n, d) < 0.0)
		{
		  faces[nbfaces].n.x *= -1;
		  faces[nbfaces].n.y *= -1;
		  faces[nbfaces].n.z *= -1;
		}

	      faces[nbfaces].element = element;
	      faces[nbfaces].physreg = -1;
	      faces[nbfaces].elemreg = -1;
              faces[nbfaces].partition = -1;
	      faces[nbfaces].pair = -1;

	      nbfaces++;

	    }

	  for (j = 2; j < 5; j++)
	    {

	      elements[element].face[j] = nbfaces;

	      faces[nbfaces].type = QUADRANGLE;
	      faces[nbfaces].nbnodes = 4;

	      switch (j)
		{
		case 2:

		  // Allocate memory
		  faces[nbfaces].node =
		    calloc (faces[nbfaces].nbnodes, sizeof (int));

		  faces[nbfaces].node[0] = elements[element].node[0];
		  faces[nbfaces].node[1] = elements[element].node[2];
		  faces[nbfaces].node[2] = elements[element].node[5];
		  faces[nbfaces].node[3] = elements[element].node[3];

		  break;

		case 3:

		  // Allocate memory
		  faces[nbfaces].node =
		    calloc (faces[nbfaces].nbnodes, sizeof (int));

		  faces[nbfaces].node[0] = elements[element].node[1];
		  faces[nbfaces].node[1] = elements[element].node[4];
		  faces[nbfaces].node[2] = elements[element].node[5];
		  faces[nbfaces].node[3] = elements[element].node[2];

		  break;

		case 4:

		  // Allocate memory
		  faces[nbfaces].node =
		    calloc (faces[nbfaces].nbnodes, sizeof (int));

		  faces[nbfaces].node[0] = elements[element].node[0];
		  faces[nbfaces].node[1] = elements[element].node[3];
		  faces[nbfaces].node[2] = elements[element].node[4];
		  faces[nbfaces].node[3] = elements[element].node[1];

		  break;

		}

	      // Calculate area 
	      faces[nbfaces].Aj =
		GeoCalcQuadArea (nodes[faces[nbfaces].node[0]],
				 nodes[faces[nbfaces].node[1]],
				 nodes[faces[nbfaces].node[2]],
				 nodes[faces[nbfaces].node[3]]);

	      // Calculate centroid   
	      faces[nbfaces].cface =
		GeoCalcCentroid4 (nodes[faces[nbfaces].node[0]],
				  nodes[faces[nbfaces].node[1]],
				  nodes[faces[nbfaces].node[2]],
				  nodes[faces[nbfaces].node[3]]);

	      // Calculate normal of the face                         
	      faces[nbfaces].n = GeoCalcNormal (nodes[faces[nbfaces].node[0]],
						nodes[faces[nbfaces].node[1]],
						nodes[faces[nbfaces].
						      node[2]]);

	      d =
		GeoSubVectorVector (faces[nbfaces].cface,
				    elements[element].celement);

	      // Normal should point to neighbour
	      // Flip normal if necessary
	      if (GeoDotVectorVector (faces[nbfaces].n, d) < 0.0)
		{
		  faces[nbfaces].n.x *= -1;
		  faces[nbfaces].n.y *= -1;
		  faces[nbfaces].n.z *= -1;
		}

	      faces[nbfaces].element = element;
	      faces[nbfaces].physreg = -1;
	      faces[nbfaces].elemreg = -1;
              faces[nbfaces].partition = -1;
	      faces[nbfaces].pair = -1;
	      faces[nbfaces].bc = NONE;

	      nbfaces++;

	    }

	}

    }

  faces = realloc (faces, nbfaces * sizeof (msh_face));

}

void
MshAddBoundaryFaces ()
{

  int i;

  nbfaces += nbpatches;

  faces = realloc (faces, nbfaces * sizeof (msh_face));

  for (i = nbfaces - nbpatches; i < nbfaces; i++)
    {

      faces[i] = patches[i - nbfaces + nbpatches];

    }

}

void
MshConnectFaces ()
{

  int i, j, k;

  double min[3], max[3];

  int nb;
  int nbv1, nbv2;
  int face, another_face;
  double dist;
  msh_vector cent[2];

  oct_data *tab = NULL;

  // Find pairs of all faces 
  // Boundary faces have no pair: -1
  // Create an octree with face cfaces

  if (nbfaces == 0)
    return;

  tab = malloc (nbfaces * sizeof (oct_data));

  min[0] = faces[0].cface.x;
  min[1] = faces[0].cface.y;
  min[2] = faces[0].cface.z;

  max[0] = faces[0].cface.x;
  max[1] = faces[0].cface.y;
  max[2] = faces[0].cface.z;

  for (i = 0; i < nbfaces; i++)
    {

      tab[i].x = faces[i].cface.x;
      tab[i].y = faces[i].cface.y;
      tab[i].z = faces[i].cface.z;

      min[0] = LMIN (min[0], faces[i].cface.x);
      min[1] = LMIN (min[1], faces[i].cface.y);
      min[2] = LMIN (min[2], faces[i].cface.z);

      max[0] = LMAX (max[0], faces[i].cface.x);
      max[1] = LMAX (max[1], faces[i].cface.y);
      max[2] = LMAX (max[2], faces[i].cface.z);

    }

  OctCreateOctree (min, max, tab, nbfaces);

  for (i = 0; i < nbleafs; i++)
    {

      nb = leafs[i].nbentities;

      // For each node, we find the two corresponding faces
      for (j = 0; j < nb; j++)
	{

	  face = leafs[i].entities[j];
	  nbv1 = faces[face].nbnodes;
	  cent[0] = faces[face].cface;

	  if (faces[face].pair != -1)
	    continue;

	  for (k = j + 1; k < nb; k++)
	    {

	      another_face = leafs[i].entities[k];
	      nbv2 = faces[another_face].nbnodes;
	      cent[1] = faces[another_face].cface;

	      dist = GeoMagVector (GeoSubVectorVector (cent[0], cent[1]));

	      if (dist < EPSILON && nbv1 == nbv2 && face != another_face)
		{

		  faces[face].pair = another_face;
		  faces[another_face].pair = face;

		  break;
		}


	    }

	  leafs[i].entities[k] = leafs[i].entities[nb - j - 1];

	}
    }

  OctDestroyOctree ();

  free (tab);

}

void
MshRemoveBoundaryFaces ()
{

  int i;

  int face, pair;

  nbfaces -= nbpatches;

  for (i = 0; i < nbfaces; i++)
    {
      face = i;

      pair = faces[face].pair;

      if (pair != -1)
	{
	  faces[face].physreg = faces[pair].physreg;
	  faces[face].elemreg = faces[pair].elemreg;
	  faces[face].partition = faces[pair].partition;
	  faces[face].bc = faces[pair].bc;
	}
    }

  for (i = 0; i < nbfaces; i++)
    {
      face = i;

      faces[face].pair = -1;
    }

  faces = realloc (faces, nbfaces * sizeof (msh_face));

}

void
MshCalcNewElementCenters ()
{


  int i, j;


  int element, face;


  double x[8], y[8], z[8];


  for (i = 0; i < nbelements; i++)

    {


      element = i;


      if (elements[element].type == TETRAHEDRON)

	{


	  // Element centre = centre of sphere connecting face centers

	  elements[element].celement.x = 0.0;

	  elements[element].celement.y = 0.0;

	  elements[element].celement.z = 0.0;


	  for (j = 0; j < elements[element].nbfaces; j++)

	    {


	      face = elements[element].face[j];


	      x[j] = faces[face].cface.x;

	      y[j] = faces[face].cface.y;

	      z[j] = faces[face].cface.z;


	    }


	  elements[element].celement.x =
	    -0.5 * (y[2] * z[3] - y[2] * z[0] - y[0] * z[3] - z[2] * y[3] +
		    z[2] * y[0] + z[0] * y[3]) / (-x[2] * y[1] * z[0] +
						  x[1] * y[0] * z[3] -
						  x[1] * z[0] * y[3] +
						  x[2] * z[1] * y[0] -
						  x[3] * y[1] * z[2] -
						  x[1] * y[2] * z[3] +
						  x[1] * y[2] * z[0] +
						  x[1] * z[2] * y[3] -
						  x[1] * z[2] * y[0] +
						  x[0] * y[2] * z[3] -
						  x[0] * z[2] * y[3] +
						  x[2] * y[1] * z[3] -
						  x[2] * y[0] * z[3] -
						  x[2] * z[1] * y[3] +
						  x[2] * z[0] * y[3] -
						  x[0] * y[1] * z[3] +
						  x[0] * z[1] * y[3] +
						  x[3] * y[1] * z[0] +
						  x[3] * z[2] * y[0] +
						  x[3] * z[1] * y[2] -
						  x[3] * z[1] * y[0] -
						  x[3] * y[2] * z[0] +
						  x[0] * y[1] * z[2] -
						  x[0] * z[1] * y[2]) *
	    (x[1] * x[1] - x[0] * x[0] + y[1] * y[1] - y[0] * y[0] +
	     z[1] * z[1] - z[0] * z[0]) + 0.5 * (-y[1] * z[0] + z[1] * y[0] +
						 y[1] * z[3] - y[0] * z[3] -
						 z[1] * y[3] +
						 z[0] * y[3]) / (-x[2] *
								 y[1] * z[0] +
								 x[1] * y[0] *
								 z[3] -
								 x[1] * z[0] *
								 y[3] +
								 x[2] * z[1] *
								 y[0] -
								 x[3] * y[1] *
								 z[2] -
								 x[1] * y[2] *
								 z[3] +
								 x[1] * y[2] *
								 z[0] +
								 x[1] * z[2] *
								 y[3] -
								 x[1] * z[2] *
								 y[0] +
								 x[0] * y[2] *
								 z[3] -
								 x[0] * z[2] *
								 y[3] +
								 x[2] * y[1] *
								 z[3] -
								 x[2] * y[0] *
								 z[3] -
								 x[2] * z[1] *
								 y[3] +
								 x[2] * z[0] *
								 y[3] -
								 x[0] * y[1] *
								 z[3] +
								 x[0] * z[1] *
								 y[3] +
								 x[3] * y[1] *
								 z[0] +
								 x[3] * z[2] *
								 y[0] +
								 x[3] * z[1] *
								 y[2] -
								 x[3] * z[1] *
								 y[0] -
								 x[3] * y[2] *
								 z[0] +
								 x[0] * y[1] *
								 z[2] -
								 x[0] * z[1] *
								 y[2]) *
	    (x[2] * x[2] - x[0] * x[0] + y[2] * y[2] - y[0] * y[0] +
	     z[2] * z[2] - z[0] * z[0]) - 0.5 * (+y[1] * z[2] - y[1] * z[0] -
						 z[2] * y[0] - z[1] * y[2] +
						 z[1] * y[0] +
						 y[2] * z[0]) / (-x[2] *
								 y[1] * z[0] +
								 x[1] * y[0] *
								 z[3] -
								 x[1] * z[0] *
								 y[3] +
								 x[2] * z[1] *
								 y[0] -
								 x[3] * y[1] *
								 z[2] -
								 x[1] * y[2] *
								 z[3] +
								 x[1] * y[2] *
								 z[0] +
								 x[1] * z[2] *
								 y[3] -
								 x[1] * z[2] *
								 y[0] +
								 x[0] * y[2] *
								 z[3] -
								 x[0] * z[2] *
								 y[3] +
								 x[2] * y[1] *
								 z[3] -
								 x[2] * y[0] *
								 z[3] -
								 x[2] * z[1] *
								 y[3] +
								 x[2] * z[0] *
								 y[3] -
								 x[0] * y[1] *
								 z[3] +
								 x[0] * z[1] *
								 y[3] +
								 x[3] * y[1] *
								 z[0] +
								 x[3] * z[2] *
								 y[0] +
								 x[3] * z[1] *
								 y[2] -
								 x[3] * z[1] *
								 y[0] -
								 x[3] * y[2] *
								 z[0] +
								 x[0] * y[1] *
								 z[2] -
								 x[0] * z[1] *
								 y[2]) *
	    (x[3] * x[3] - x[0] * x[0] + y[3] * y[3] - y[0] * y[0] +
	     z[3] * z[3] - z[0] * z[0]);

	  elements[element].celement.y =
	    +0.5 * (x[2] * z[3] - x[2] * z[0] - x[0] * z[3] - z[2] * x[3] +
		    z[2] * x[0] + z[0] * x[3]) / (-x[2] * y[1] * z[0] +
						  x[1] * y[0] * z[3] -
						  x[1] * z[0] * y[3] +
						  x[2] * z[1] * y[0] -
						  x[3] * y[1] * z[2] -
						  x[1] * y[2] * z[3] +
						  x[1] * y[2] * z[0] +
						  x[1] * z[2] * y[3] -
						  x[1] * z[2] * y[0] +
						  x[0] * y[2] * z[3] -
						  x[0] * z[2] * y[3] +
						  x[2] * y[1] * z[3] -
						  x[2] * y[0] * z[3] -
						  x[2] * z[1] * y[3] +
						  x[2] * z[0] * y[3] -
						  x[0] * y[1] * z[3] +
						  x[0] * z[1] * y[3] +
						  x[3] * y[1] * z[0] +
						  x[3] * z[2] * y[0] +
						  x[3] * z[1] * y[2] -
						  x[3] * z[1] * y[0] -
						  x[3] * y[2] * z[0] +
						  x[0] * y[1] * z[2] -
						  x[0] * z[1] * y[2]) *
	    (x[1] * x[1] - x[0] * x[0] + y[1] * y[1] - y[0] * y[0] +
	     z[1] * z[1] - z[0] * z[0]) - 0.5 * (+x[1] * z[3] - x[1] * z[0] -
						 x[0] * z[3] - z[1] * x[3] +
						 z[1] * x[0] +
						 z[0] * x[3]) / (-x[2] *
								 y[1] * z[0] +
								 x[1] * y[0] *
								 z[3] -
								 x[1] * z[0] *
								 y[3] +
								 x[2] * z[1] *
								 y[0] -
								 x[3] * y[1] *
								 z[2] -
								 x[1] * y[2] *
								 z[3] +
								 x[1] * y[2] *
								 z[0] +
								 x[1] * z[2] *
								 y[3] -
								 x[1] * z[2] *
								 y[0] +
								 x[0] * y[2] *
								 z[3] -
								 x[0] * z[2] *
								 y[3] +
								 x[2] * y[1] *
								 z[3] -
								 x[2] * y[0] *
								 z[3] -
								 x[2] * z[1] *
								 y[3] +
								 x[2] * z[0] *
								 y[3] -
								 x[0] * y[1] *
								 z[3] +
								 x[0] * z[1] *
								 y[3] +
								 x[3] * y[1] *
								 z[0] +
								 x[3] * z[2] *
								 y[0] +
								 x[3] * z[1] *
								 y[2] -
								 x[3] * z[1] *
								 y[0] -
								 x[3] * y[2] *
								 z[0] +
								 x[0] * y[1] *
								 z[2] -
								 x[0] * z[1] *