mesh.c

OpenFVM-v1.1 open source cfd code

								 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 * (-x[1] * z[0] + x[1] * z[2] -
						 z[2] * x[0] - z[1] * x[2] +
						 x[2] * z[0] +
						 z[1] * x[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.z =
	    -0.5 * (x[2] * y[3] - x[2] * y[0] - x[0] * y[3] - y[2] * x[3] +
		    y[2] * x[0] + y[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 * (-y[1] * x[3] + x[1] * y[3] -
						 x[1] * y[0] - x[0] * y[3] +
						 y[0] * x[3] +
						 y[1] * x[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[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 * (+x[1] * y[2] - x[1] * y[0] -
						 y[2] * x[0] - y[1] * x[2] +
						 y[1] * x[0] +
						 x[2] * y[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]);


	}


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

	{


	  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 += x[j];

	      elements[element].celement.y += y[j];

	      elements[element].celement.z += z[j];


	    }


	  elements[element].celement.x /= elements[element].nbfaces;

	  elements[element].celement.y /= elements[element].nbfaces;

	  elements[element].celement.z /= elements[element].nbfaces;


	}


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

	{


	  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 += x[j];

	      elements[element].celement.y += y[j];

	      elements[element].celement.z += z[j];


	    }


	  elements[element].celement.x /= elements[element].nbfaces;

	  elements[element].celement.y /= elements[element].nbfaces;

	  elements[element].celement.z /= elements[element].nbfaces;


	}


    }


}

void
MshCalcPropMesh ()
{

  int i, j;

  int element, face, pair;

  double sum_area;
  double sum_volume;

  double total_area;
  double total_volume;

  total_area = 0.0;
  total_volume = 0.0;

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

      element = i;

      sum_area = 0.0;
      sum_volume = 0.0;

      // Volume calculation using Gauss theorem

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

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

	  sum_volume += faces[face].cface.x * faces[face].A.x;

	  pair = faces[face].pair;

	  if (pair == -1)
	    sum_area += faces[face].Aj;

	}

      if (sum_volume <= 0.0)
	{
	  PetscPrintf (PETSC_COMM_WORLD,
		       "\nError: Element with zero or negative volume\n");
	  exit (LOGICAL_ERROR);
	}

      elements[element].Vp = sum_volume;

      total_volume += elements[i].Vp;
      total_area += sum_area;

    }

  PetscPrintf (PETSC_COMM_WORLD, "Sub total surface area: \t\t\t%+.3E %s^2\n",
	       total_area, parameter.ulength);
  PetscPrintf (PETSC_COMM_WORLD, "Sub total volume: \t\t\t\t%+.3E %s^3\n",
	       total_volume, parameter.ulength);

}

void
MshCalcPropFaces ()
{

  int i, j, k;

  int element, face, pair, neighbor;

  int nb_internal_faces;
  int nb_total_faces;

  double distance_mean;
  double northo_mean, shapedev_mean;

  int node1, node2;
  double d, dmin, dmax;
  double aspectratio, aspectratio_min, aspectratio_max;

  nb_internal_faces = 0;
  nb_total_faces = 0;

  northo_mean = 0.0;
  distance_mean = 0.0;
  shapedev_mean = 0.0;

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

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      faces[face].rpl = GeoSubVectorVector (faces[face].cface, GeoMultScalarVector (GeoDotVectorVector(GeoSubVectorVector(faces[face].cface, elements[element].celement), faces[face].n), faces[face].n));

      if (pair != -1)
	{

	  neighbor = faces[pair].element;

	  faces[face].rnl =
	    GeoSubVectorVector (faces[face].cface, GeoMultScalarVector (GeoDotVectorVector(GeoSubVectorVector(faces[face].cface, elements[neighbor].celement), faces[face].n), faces[face].n));

	  faces[face].A = GeoMultScalarVector (faces[face].Aj, faces[face].n);

	  faces[face].d = GeoSubVectorVector (faces[face].rnl, faces[face].rpl);

	  faces[face].dj = GeoMagVector (faces[face].d);

          if (faces[face].dj == 0)
	  {
		  PetscPrintf (PETSC_COMM_WORLD, "\nError: Problem with mesh\n");
		  exit (LOGICAL_ERROR);
	  }

	  distance_mean += faces[face].dj;
	  nb_internal_faces++;

	}
      else
	{

	  faces[face].A = GeoMultScalarVector (faces[face].Aj, faces[face].n);

	  faces[face].d = GeoSubVectorVector (faces[face].cface, faces[face].rpl);

	  faces[face].dj = GeoMagVector (faces[face].d);

          if (faces[face].dj == 0)
	  {
		  PetscPrintf (PETSC_COMM_WORLD, "\nError: Problem with mesh\n");
		  exit (LOGICAL_ERROR);
	  }

	}

      // Measure non-orthogonality of the mesh
      northo_mean +=
	GeoMagVector (GeoSubVectorVector
		      (elements[element].celement,
		       faces[face].rpl)) / faces[face].dj;
      nb_total_faces++;

    }

  if (nb_internal_faces > 0)
    distance_mean /= nb_internal_faces;

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

      pair = faces[face].pair;

      if (pair != -1)
	{
	  shapedev_mean +=
	    LABS (faces[face].dj - distance_mean) / distance_mean;

	}

    }

  if (nb_total_faces > 0)
    northo_mean /= nb_total_faces;

  if (nb_internal_faces > 0)
    shapedev_mean /= nb_internal_faces;

  aspectratio_min = GREAT;
  aspectratio_max = SMALL;

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

      element = i;

      dmin = 1E+30f;
      dmax = 1E-30f;

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

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

	  for (k = 0; k < faces[face].nbnodes; k++)
	    {

	      node1 = faces[face].node[(k + 0) % faces[face].nbnodes];
	      node2 = faces[face].node[(k + 1) % faces[face].nbnodes];

	      d =
		(float)
		GeoMagVector (GeoSubVectorVector
			      (nodes[node2], nodes[node1]));

	      dmin = LMIN (dmin, d);
	      dmax = LMAX (dmax, d);

	    }

	}

      aspectratio = dmin / dmax;

      aspectratio_min = LMIN (aspectratio_min, aspectratio);
      aspectratio_max = LMAX (aspectratio_max, aspectratio);

    }

  PetscPrintf (PETSC_COMM_WORLD, "\n");
  PetscPrintf (PETSC_COMM_WORLD, "Non-orthogonality of the mesh: %f\n",
	       northo_mean);
  PetscPrintf (PETSC_COMM_WORLD, "Shape deviation: %f\n", shapedev_mean);
  PetscPrintf (PETSC_COMM_WORLD, "Aspect ratio (worst - best): %f - %f\n",
	       aspectratio_min, aspectratio_max);

}

int
MshImportMSH (char *file)
{
  int i, j;

  int nnod;
  int nele;

  int eindex;
  int maxindex;
  int etype;

  int ntags;
  int physreg, elemreg, partition;

  int ival;

  float fval;

  FILE *fp;
  char descr[512];

  fp = fopen (file, "r");

  if (fp == NULL)
    {
      PetscPrintf (PETSC_COMM_WORLD, "\nError: Mesh file not found!\n");
      PetscPrintf (PETSC_COMM_WORLD, "%s\n\n", file);
      exit (LOGICAL_ERROR);
    }

  PetscPrintf (PETSC_COMM_WORLD, "\nReading mesh file: %s ...\n", file);

  maxindex = 0;

  do
    {

      do
	{

	  strcpy (descr, "");

	  fscanf (fp, "%s", descr);

	  if (strcmp (descr, "$Nodes") == 0)
	    break;

	  if (strcmp (descr, "$Elements") == 0)
	    break;

	}
      while (!feof (fp));

      if (strcmp (descr, "$Nodes") == 0)
	{
	  fscanf (fp, "%d", &nnod);

	  nodes = realloc (nodes, nnod * sizeof (msh_vector));

	  nod_correlation_malloced = nnod;
	  nod_correlation =
	    realloc (nod_correlation,
		     nod_correlation_malloced * sizeof (int));

	  nbnodes = 0;

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

	      fscanf (fp, "%d", &ival);

	      if (ival >= nod_correlation_malloced)
		{
		  nod_correlation_malloced *= 2;
		  nod_correlation =
		    realloc (nod_correlation,
			     nod_correlation_malloced * sizeof (int));
		}

	      nod_correlation[ival] = i;

	      fscanf (fp, "%f", &fval);
	      nodes[nbnodes].x = fval;

	      fscanf (fp, "%f", &fval);
	      nodes[nbnodes].y = fval;

	      fscanf (fp, "%f", &fval);
	      nodes[nbnodes].z = fval;

	      nbnodes++;

	    }
	}

      if (strcmp (descr, "$Elements") == 0)
	{

	  fscanf (fp, "%d", &nele);

	  patches = realloc (patches, nele * sizeof (msh_face));
	  elements = realloc (elements, nele * sizeof (msh_element));

	  nbpatches = 0;
	  nbelements = 0;

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

	      fscanf (fp, "%d %d", &eindex, &etype);

	      if (etype < 1)
		{
		  // Invalid element

		  PetscPrintf (PETSC_COMM_WORLD,
			       "\nError: Invalid element\n");
		  exit (LOGICAL_ERROR);

		}

	      if (etype == 1)
		{

		  // Beam element

		  PetscPrintf (PETSC_COMM_WORLD,
			       "\nError: Invalid element\n");
		  exit (LOGICAL_ERROR);

		}

	      if (etype == 2 || etype == 3)