mesh.c

OpenFVM-v1.1 open source cfd code

/***************************************************************************
 *   Copyright (C) 2004-2008 by OpenFVM team                               *
 *   http://sourceforge.net/projects/openfvm/                              *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>

// PETSc
#include "petscksp.h"

#include "globals.h"
#include "param.h"
#include "mesh.h"
#include "bcond.h"
#include "geocalc.h"
#include "octree.h"

void
MshFreeMemory ()
{

  int i;

  int patch, face, element;

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

      face = i;

      if (faces[face].type == TRIANGLE || faces[face].type == QUADRANGLE)
	{
	  free (faces[face].node);
	}

    }

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

      free (patches[patch].node);

    }

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

      element = i;

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

	  free (elements[element].b);
	  free (elements[element].c);
	  free (elements[element].d);

	}

    }

  if (nbnodes > 0)
    {
      nbnodes = 0;
      free (nodes);
    }

  if (nbfaces > 0)
    {
      nbfaces = 0;
      free (faces);
    }

  if (nbelements > 0)
    {
      nbelements = 0;
      free (elements);
    }

  if (nbpatches > 0)
    {
      nbpatches = 0;
      free (patches);
    }

}

void
MshCorrectNonOrthogonality ()
{

  int i;

  int face, pair;

  int element;
  
  //int neighbor;

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

      face = i;

      element = faces[face].element;

      pair = faces[face].pair;

      faces[face].kj = 0.0;

      /*
      if (pair != -1)
	{

	  neighbor = faces[pair].element;

	  faces[face].kj +=
	    GeoMagVector (GeoSubVectorVector
			  (elements[element].celement, faces[face].rpl));

	  faces[face].kj +=
	    GeoMagVector (GeoSubVectorVector
			  (elements[neighbor].celement, faces[pair].rpl));


	}
      */
    }


}

void
MshGetShapeFunctions ()
{

  int i;

  int element;

  //double detJ;

  double x14, y14, z14;
  double x24, y24, z24;
  double x34, y34, z34;

  double bi, bj, bk, bl;
  double ci, cj, ck, cl;
  double di, dj, dk, dl;

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

      element = i;

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

	  x14 =
	    nodes[elements[element].node[0]].x -
	    nodes[elements[element].node[3]].x;
	  y14 =
	    nodes[elements[element].node[0]].y -
	    nodes[elements[element].node[3]].y;
	  z14 =
	    nodes[elements[element].node[0]].z -
	    nodes[elements[element].node[3]].z;

	  x24 =
	    nodes[elements[element].node[1]].x -
	    nodes[elements[element].node[3]].x;
	  y24 =
	    nodes[elements[element].node[1]].y -
	    nodes[elements[element].node[3]].y;
	  z24 =
	    nodes[elements[element].node[1]].z -
	    nodes[elements[element].node[3]].z;

	  x34 =
	    nodes[elements[element].node[2]].x -
	    nodes[elements[element].node[3]].x;
	  y34 =
	    nodes[elements[element].node[2]].y -
	    nodes[elements[element].node[3]].y;
	  z34 =
	    nodes[elements[element].node[2]].z -
	    nodes[elements[element].node[3]].z;

	  /*
	     detJ = (x14 * y24 * z34 + y14 * z24 * x34 + z14 * x24 * y34) - 
	     (x34 * y24 * z14 + y34 * z24 * x14 + z34 * x24 * y14);

	     elements[element].Vp = LABS(detJ) / 6.0;
	   */

	  bi = (y24 * z34 - y34 * z24);
	  bj = (y34 * z14 - y14 * z34);
	  bk = (y14 * z24 - y24 * z14);
	  bl = -(bi + bj + bk);

	  ci = (z24 * x34 - z34 * x24);
	  cj = (z34 * x14 - z14 * x34);
	  ck = (z14 * x24 - z24 * x14);
	  cl = -(ci + cj + ck);

	  di = (x24 * y34 - x34 * y24);
	  dj = (x34 * y14 - x14 * y34);
	  dk = (x14 * y24 - x24 * y14);
	  dl = -(di + dj + dk);

	  // Allocate memory

	  elements[element].b =
	    calloc (elements[element].nbnodes, sizeof (double));
	  elements[element].c =
	    calloc (elements[element].nbnodes, sizeof (double));
	  elements[element].d =
	    calloc (elements[element].nbnodes, sizeof (double));

	  elements[element].b[0] = bi;
	  elements[element].b[1] = bj;
	  elements[element].b[2] = bk;
	  elements[element].b[3] = bl;

	  elements[element].c[0] = ci;
	  elements[element].c[1] = cj;
	  elements[element].c[2] = ck;
	  elements[element].c[3] = cl;

	  elements[element].d[0] = di;
	  elements[element].d[1] = dj;
	  elements[element].d[2] = dk;
	  elements[element].d[3] = dl;

	}
    }

}

void
MshGetElementTypes ()
{

  int i;

  int element, face, pair;

  nbtris = 0;
  nbquads = 0;
  nbtetras = 0;
  nbhexas = 0;
  nbprisms = 0;

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

      face = i;

      pair = faces[face].pair;

      if (faces[face].type == TRIANGLE)
	{
	  nbtris++;
	}

      if (faces[face].type == QUADRANGLE)
	{
	  nbquads++;
	}

    }

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

      element = i;

      if (elements[element].type == TETRAHEDRON)
	{
	  nbtetras++;
	}

      if (elements[element].type == HEXAHEDRON)
	{
	  nbhexas++;
	}

      if (elements[element].type == PRISM)
	{
	  nbprisms++;
	}
    }

}

void
MshCreateFaces ()
{

  int i, j;

  int element;

  msh_vector d;

  // Create faces
  faces = realloc (faces, nbelements * 8 * sizeof (msh_face));

  nbfaces = 0;

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

      element = i;

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

	  elements[element].nbfaces = 1;

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

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

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

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

	  // 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]]);

	  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 == QUADRANGLE)
	{

	  elements[element].nbfaces = 1;

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

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

	  faces[nbfaces].type = QUADRANGLE;
	  faces[nbfaces].nbnodes = 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[1];
	  faces[nbfaces].node[2] = elements[element].node[2];
	  faces[nbfaces].node[3] = elements[element].node[3];

	  // 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]]);

	  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 == TETRAHEDRON)
	{

	  elements[element].nbfaces = 4;

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

	  for (j = 0; j < 4; 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[1];
		  faces[nbfaces].node[1] = elements[element].node[3];
		  faces[nbfaces].node[2] = elements[element].node[2];

		  break;

		case 2:

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

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

		  break;

		case 3:

		  // 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[1];
		  faces[nbfaces].node[2] = elements[element].node[0];

		  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;
	      faces[nbfaces].bc = NONE;

	      nbfaces++;

	    }

	}

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

	  elements[element].nbfaces = 6;

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

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

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

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

	      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];
		  faces[nbfaces].node[3] = elements[element].node[3];

		  break;

		case 1:

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

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

		  break;

		case 2:

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

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

		  break;

		case 3:

		  // 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[7];
		  faces[nbfaces].node[2] = elements[element].node[4];
		  faces[nbfaces].node[3] = elements[element].node[0];

		  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[4];
		  faces[nbfaces].node[2] = elements[element].node[5];
		  faces[nbfaces].node[3] = elements[element].node[1];

		  break;

		case 5:

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

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

		  break;

		}

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