mesh.c

OpenFVM-v1.1 open source cfd code

		{

		  // Shell element or surface element

		  if (etype == 2) ival = 3;
		  if (etype == 3) ival = 4;

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

		  if (ntags != 3)
		  {
			printf ("\nError: Number of tags != 3\n");
		 	exit (LOGICAL_ERROR);
		  }

		  fscanf (fp, "%d %d %d", &physreg, &elemreg, &partition);

		  patches[nbpatches].nbnodes = ival;

		  patches[nbpatches].cface.x = 0.0;
		  patches[nbpatches].cface.y = 0.0;
		  patches[nbpatches].cface.z = 0.0;

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

		  for (j = 0; j < patches[nbpatches].nbnodes; j++)
		    {
		      fscanf (fp, "%d", &ival);

		      patches[nbpatches].node[j] = nod_correlation[ival];

		      patches[nbpatches].cface.x +=
			nodes[patches[nbpatches].node[j]].x;
		      patches[nbpatches].cface.y +=
			nodes[patches[nbpatches].node[j]].y;
		      patches[nbpatches].cface.z +=
			nodes[patches[nbpatches].node[j]].z;

		    }

		  patches[nbpatches].cface.x /= patches[nbpatches].nbnodes;
		  patches[nbpatches].cface.y /= patches[nbpatches].nbnodes;
		  patches[nbpatches].cface.z /= patches[nbpatches].nbnodes;

		  patches[nbpatches].physreg = physreg;
		  patches[nbpatches].elemreg = elemreg;
		  patches[nbpatches].partition = partition;

		  patches[nbpatches].element = -1;
		  patches[nbpatches].pair = -1;
		  patches[nbpatches].bc = NONE;

		  if (etype == 2)
		    patches[nbpatches].type = TRIANGLE;

		  if (etype == 3)
		    patches[nbpatches].type = QUADRANGLE;

		  patches[nbpatches].index = eindex;

		  if (nbelements == 0)
		    maxindex = eindex;

		  nbpatches++;

		}

	      if (etype >= 4 && etype <= 6)
		{

		  // Solid elements or volume element

		  if (etype == 4) ival = 4;
		  if (etype == 5) ival = 8;
		  if (etype == 6) ival = 6;

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

		  if (ntags != 3)
		  {
			printf ("\nError: Number of tags != 3\n");
		 	exit (LOGICAL_ERROR);
		  }

		  fscanf (fp, "%d %d %d", &physreg, &elemreg, &partition);

		  elements[nbelements].nbnodes = ival;

		  elements[nbelements].celement.x = 0.0;
		  elements[nbelements].celement.y = 0.0;
		  elements[nbelements].celement.z = 0.0;

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

		  for (j = 0; j < elements[nbelements].nbnodes; j++)
		    {
		      fscanf (fp, "%d", &ival);

		      elements[nbelements].node[j] = nod_correlation[ival];

		      elements[nbelements].celement.x +=
			nodes[elements[nbelements].node[j]].x;
		      elements[nbelements].celement.y +=
			nodes[elements[nbelements].node[j]].y;
		      elements[nbelements].celement.z +=
			nodes[elements[nbelements].node[j]].z;

		    }

		  elements[nbelements].celement.x /= elements[nbelements].
		    nbnodes;
		  elements[nbelements].celement.y /= elements[nbelements].
		    nbnodes;
		  elements[nbelements].celement.z /= elements[nbelements].
		    nbnodes;

		  elements[nbelements].physreg = physreg;
		  elements[nbelements].elemreg = elemreg;
		  elements[nbelements].partition = partition;

		  if (etype == 4)
		    elements[nbelements].type = TETRAHEDRON;

		  if (etype == 5)
		    elements[nbelements].type = HEXAHEDRON;

		  if (etype == 6)
		    elements[nbelements].type = PRISM;

		  elements[nbelements].index = eindex;

		  nbelements++;

		}

	      if (etype >= 7)
		{

		  // Unkown element

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

		}

	    }

	}

    }
  while (!feof (fp));

  fclose (fp);

  nodes = realloc (nodes, nbnodes * sizeof (msh_vector));
  patches = realloc (patches, nbpatches * sizeof (msh_face));
  elements = realloc (elements, nbelements * sizeof (msh_element));

  free (nod_correlation);

  // Assign patches for parallel processing

  for (i = 0; i < nbpatches; i++)
    {
      if (patches[i].index > maxindex)
	{
	  patches[i].bc = PROCESSOR;
	}
    }

  PetscPrintf (PETSC_COMM_WORLD, "Done.\n");

  MshCreateFaces ();
  MshAddBoundaryFaces ();
  MshConnectFaces ();
  MshRemoveBoundaryFaces ();
  MshConnectFaces ();
  MshCalcPropFaces ();
  //MshCalcNewElementCenters();

  PetscPrintf (PETSC_COMM_WORLD, "\n");
  PetscPrintf (PETSC_COMM_WORLD, "Number of nodes: \t\t\t\t%d\n", nbnodes);
  PetscPrintf (PETSC_COMM_WORLD, "Number of faces: \t\t\t\t%d\n", nbfaces);
  PetscPrintf (PETSC_COMM_WORLD, "Number of boundary faces: \t\t\t%d\n",
	       nbpatches);
  PetscPrintf (PETSC_COMM_WORLD, "Number of elements: \t\t\t\t%d\n",
	       nbelements);

  if (nbelements == 0)
    {
      PetscPrintf (PETSC_COMM_WORLD, "\nError: No solid elements.\n");
      exit (LOGICAL_ERROR);
    }

  MshCalcPropMesh ();
  MshGetElementTypes ();
  MshCorrectNonOrthogonality ();
  //MshGetShapeFunctions();

  return LOGICAL_TRUE;

}

int
MshExportMSH (char *file)
{

  int i, j;

  int node, patch, element;

  FILE *fp;

  fp = fopen (file, "w");

  if (fp == NULL)
    return LOGICAL_FALSE;

  fprintf (fp, "$MeshFormat\n");
  fprintf (fp, "2 0 8\n");
  fprintf (fp, "$EndMeshFormat\n");

  fprintf (fp, "$Nodes\n");
  fprintf (fp, "%d\n", nbnodes);

  for (i = 0; i < nbnodes; i++)
    {
      node = i;

      fprintf (fp, "%d %f %f %f", node + 1, nodes[node].x, nodes[node].y, nodes[node].z);
      fprintf (fp, "\n");
    }

  fprintf (fp, "$EndNodes\n");

  fprintf (fp, "$Elements\n");

  fprintf (fp, "%d\n", nbpatches + nbelements);

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

      patch = i;

      fprintf (fp, "%d", patches[patch].index);

      if (patches[patch].type == TRIANGLE)
	fprintf (fp, " %d", 2);

      if (patches[patch].type == QUADRANGLE)
	fprintf (fp, " %d", 3);

      fprintf (fp, " %d", 3);

      fprintf (fp, " %d", patches[patch].physreg);
      fprintf (fp, " %d", patches[patch].elemreg);
      fprintf (fp, " %d", patches[patch].partition);

      for (j = 0; j < patches[patch].nbnodes; j++)
	{
	  fprintf (fp, " %d", patches[patch].node[j] + 1);
	}

      fprintf (fp, "\n");
    }

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

      element = i;

      fprintf (fp, "%d", elements[element].index);

      if (elements[element].type == TETRAHEDRON)
	fprintf (fp, " %d", 4);

      if (elements[element].type == HEXAHEDRON)
	fprintf (fp, " %d", 5);

      if (elements[element].type == PRISM)
	fprintf (fp, " %d", 6);

      fprintf (fp, " %d", 3);

      fprintf (fp, " %d", elements[element].physreg);
      fprintf (fp, " %d", elements[element].elemreg);
      fprintf (fp, " %d", elements[element].partition);

      for (j = 0; j < elements[element].nbnodes; j++)
	{
	  fprintf (fp, " %d", elements[element].node[j] + 1);
	}

      fprintf (fp, "\n");
    }

  fprintf (fp, "$EndElements\n");

  fclose (fp);

  return LOGICAL_TRUE;

}

int
MshExportDecomposedMSH (char *file, int partition, int nbpartitions)
{

  int i, j, n;

  int node, face, pair, patch, element, neighbor;

  float b;

  FILE *fp;

  int *nflag;

  int nbnewnodes;
  int nbnewpatches;
  int nbnewelements;

  msh_vector *newnodes;
  msh_face *newpatches;
  msh_element *newelements;

  int *newindex;

  // Identify elements in partition

  newelements = calloc (nbelements, sizeof (msh_element));

  nbnewelements = 0;

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

      element = i;

      if (elements[element].partition != partition)
	continue;

      newelements[nbnewelements] = elements[element];
      nbnewelements++;

    }

  if (nbnewelements == 0)
    {
      PetscPrintf (PETSC_COMM_WORLD, "\nError: Number of elements is zero\n");
      return LOGICAL_ERROR;
    }

  b = (float) nbnewelements / (float) nbelements *nbpartitions;

  printf ("Partition: %d, Balance: %.2f\n", partition, b);

  newelements = realloc (newelements, nbnewelements * sizeof (msh_element));

  // Remove unconnected and merged nodes

  nflag = calloc (nbnodes, sizeof (int));

  for (i = 0; i < nbnodes; i++)
    {
      node = i;

      nflag[node] = LOGICAL_FALSE;
    }

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

      element = i;

      for (j = 0; j < newelements[element].nbnodes; j++)
	{
	  node = newelements[element].node[j];

	  nflag[node] = LOGICAL_TRUE;
	}

    }

  newnodes = calloc (nbnodes, sizeof (msh_vector));

  newindex = calloc (nbnodes, sizeof (int));

  nbnewnodes = 0;

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

      node = i;

      if (nflag[node] == LOGICAL_TRUE)
	{

	  newnodes[nbnewnodes] = nodes[node];

	  newindex[node] = nbnewnodes;

	  nbnewnodes++;

	}
    }

  newnodes = realloc (newnodes, nbnewnodes * sizeof (msh_vector));

  // Identify patches in partition

  newpatches = calloc (nbpatches + nbfaces, sizeof (msh_face));

  nbnewpatches = 0;

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

      patch = i;

      n = LOGICAL_TRUE;

      for (j = 0; j < patches[patch].nbnodes; j++)
	{
	  node = patches[patch].node[j];

	  if (nflag[node] == LOGICAL_FALSE)
	    n = LOGICAL_FALSE;
	}

      if (n == LOGICAL_FALSE)
	continue;

      newpatches[nbnewpatches] = patches[patch];
      nbnewpatches++;
    }

  // Add patches which belong to cut 

  n = 0;

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

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      if (pair != -1)
	{

	  neighbor = faces[pair].element;

	  if (elements[element].partition != partition)
	    continue;

	  if (elements[element].partition == elements[neighbor].partition)
	    continue;

	  newpatches[nbnewpatches] = faces[face];

	  // Set physical region pointing to neighbor
	  newpatches[nbnewpatches].physreg = elements[neighbor].index;

	  // Set element partition pointing to neighbor region
	  newpatches[nbnewpatches].partition = elements[neighbor].partition;

	  // Set index to the end of the list of elements
	  newpatches[nbnewpatches].index = nbelements + nbpatches + n + 1;

	  n++;

	  nbnewpatches++;

	}

    }

  newpatches = realloc (newpatches, nbnewpatches * sizeof (msh_face));

  // Renumber patches and elements        

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

      element = i;

      for (j = 0; j < newelements[element].nbnodes; j++)
	{
	  node = newelements[element].node[j];

	  newelements[element].node[j] = newindex[node];

	}
    }

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

      patch = i;

      for (j = 0; j < newpatches[patch].nbnodes; j++)
	{
	  node = newpatches[patch].node[j];

	  newpatches[patch].node[j] = newindex[node];

	}
    }

  fp = fopen (file, "w");

  if (fp == NULL)
    return LOGICAL_FALSE;

  fprintf (fp, "$MeshFormat\n");
  fprintf (fp, "2 0 8\n");
  fprintf (fp, "$EndMeshFormat\n");

  fprintf (fp, "$Nodes\n");
  fprintf (fp, "%d\n", nbnewnodes);

  for (i = 0; i < nbnewnodes; i++)
    {
      node = i;

      fprintf (fp, "%d %f %f %f", node + 1, newnodes[node].x,
	       newnodes[node].y, newnodes[node].z);

      fprintf (fp, "\n");
    }

  fprintf (fp, "$EndNodes\n");

  fprintf (fp, "$Elements\n");

  fprintf (fp, "%d\n", nbnewpatches + nbnewelements);

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

      patch = i;

      if (newpatches[patch].index > nbelements + nbpatches)
	continue;

      fprintf (fp, "%d", newpatches[patch].index);

      if (newpatches[patch].type == TRIANGLE)
	fprintf (fp, " %d", 2);

      if (newpatches[patch].type == QUADRANGLE)
	fprintf (fp, " %d", 3);

      fprintf (fp, " %d", 3);

      fprintf (fp, " %d", newpatches[patch].physreg);
      fprintf (fp, " %d", newpatches[patch].elemreg);
      fprintf (fp, " %d", newpatches[patch].partition);

      for (j = 0; j < newpatches[patch].nbnodes; j++)
	{
	  fprintf (fp, " %d", newpatches[patch].node[j] + 1);
	}

      fprintf (fp, "\n");
    }

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

      element = i;

      fprintf (fp, "%d", newelements[element].index);

      if (newelements[element].type == TETRAHEDRON)
	fprintf (fp, " %d", 4);

      if (newelements[element].type == HEXAHEDRON)
	fprintf (fp, " %d", 5);

      if (newelements[element].type == PRISM)
	fprintf (fp, " %d", 6);

      fprintf (fp, " %d", 3);

      fprintf (fp, " %d", newelements[element].physreg);
      fprintf (fp, " %d", newelements[element].elemreg);
      fprintf (fp, " %d", newelements[element].partition);

      for (j = 0; j < newelements[element].nbnodes; j++)
	{
	  fprintf (fp, " %d", newelements[element].node[j] + 1);
	}

      fprintf (fp, "\n");
    }

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

      patch = i;

      if (newpatches[patch].index <= nbelements + nbpatches)
	continue;

      fprintf (fp, "%d", newpatches[patch].index);

      if (newpatches[patch].type == TRIANGLE)
	fprintf (fp, " %d", 2);

      if (newpatches[patch].type == QUADRANGLE)
	fprintf (fp, " %d", 3);

      fprintf (fp, " %d", 3);

      fprintf (fp, " %d", newpatches[patch].physreg);
      fprintf (fp, " %d", newpatches[patch].elemreg);
      fprintf (fp, " %d", newpatches[patch].partition);

      for (j = 0; j < newpatches[patch].nbnodes; j++)
	{
	  fprintf (fp, " %d", newpatches[patch].node[j] + 1);
	}

      fprintf (fp, "\n");
    }

  fprintf (fp, "$EndElements\n");

  fclose (fp);

  free (nflag);
  free (newnodes);
  free (newpatches);
  free (newelements);
  free (newindex);

  return LOGICAL_TRUE;

}