post.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>

#include "variables.h"
#include "ioutils.h"
#include "globals.h"
#include "mesh.h"
#include "param.h"
#include "bcond.h"
#include "geocalc.h"
#include "gradient.h"

#include "post.h"

// Function based on streamFunction.C developed by OpenFOAM
// modified and translated to C by the OpenFVM team in 08/02/2006

void
CalculateStreamFunction (double *streamFunction)
{

  int i, j;

  int face, pair, element;

  int found, finished;

  int *visitedPoint;
  int nVisited, nVisitedOld;

  int bPointFound, pointFound;

  double currentBStream, currentStream;

  msh_vector currentBStreamPoint, currentStreamPoint;
  msh_vector edgeHat;

  double vmin, vmax;

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

  nVisited = 0;
  nVisitedOld = 0;

  finished = LOGICAL_TRUE;

  do
    {

      // Find the boundary face with zero flux. set the stream function
      // to zero on that face

      found = LOGICAL_FALSE;

      // Boundary faces

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

	  face = i;

	  pair = faces[face].pair;

	  if (pair != -1)
	    continue;

	  if (faces[face].bc == EMPTY)
	    continue;

	  if (LABS (V_GetCmp (&uf, face + 1)) < SMALL)
	    {

	      // Zero flux face found
	      found = LOGICAL_TRUE;

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

		  if (visitedPoint[faces[face].node[j]] == 1)
		    {
		      found = LOGICAL_FALSE;
		      break;
		    }

		}

	      if (found == LOGICAL_TRUE)
		{

		  //printf("Zero face: %d\n", face);              

		  for (j = 0; j < faces[face].nbnodes; j++)
		    {
		      streamFunction[faces[face].node[j]] = 0.0;
		      visitedPoint[faces[face].node[j]] = 1;
		      nVisited++;
		    }

		  break;
		}


	    }

	  if (found == LOGICAL_TRUE)
	    break;

	}

      if (found == LOGICAL_FALSE)
	{

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

	      element = i;

	      found = LOGICAL_TRUE;

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

		  if (visitedPoint[elements[element].node[j]] == 1)
		    {
		      found = LOGICAL_FALSE;
		      break;
		    }

		}

	      if (found == LOGICAL_TRUE)
		{

		  for (j = 0; j < elements[element].nbnodes; j++)
		    {
		      streamFunction[elements[element].node[j]] = 0.0;
		      visitedPoint[elements[element].node[j]] = 1;
		      nVisited++;
		    }

		  break;
		}


	    }

	}

      // Loop through all faces. If one of the points on
      // the face has the streamFunction value different
      // from -1, all points with -1 on that face have the
      // streamFunction value equal to the face flux in
      // that point plus the value in the visited point

      do
	{

	  finished = LOGICAL_TRUE;

	  // Boundary faces

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

	      face = i;

	      pair = faces[face].pair;

	      if (pair != -1)
		continue;

	      bPointFound = LOGICAL_FALSE;

	      currentBStream = 0.0;

	      currentBStreamPoint.x = 0.0;
	      currentBStreamPoint.y = 0.0;
	      currentBStreamPoint.z = 0.0;

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

		  // Check if the point has been visited

		  if (visitedPoint[faces[face].node[j]] == 1)
		    {
		      // The point has been visited
		      currentBStream = streamFunction[faces[face].node[j]];
		      currentBStreamPoint = nodes[faces[face].node[j]];

		      bPointFound = LOGICAL_TRUE;

		      break;
		    }


		}

	      if (bPointFound == LOGICAL_TRUE)
		{

		  // Sort out other points on the face                            
		  for (j = 0; j < faces[face].nbnodes; j++)
		    {

		      // Check if the point has been visited
		      if (visitedPoint[faces[face].node[j]] == 0)
			{

			  if (faces[face].bc != EMPTY)
			    {

			      edgeHat =
				GeoSubVectorVector (nodes
						    [faces[face].node[j]],
						    currentBStreamPoint);

			      edgeHat.z = 0.0;

			      edgeHat = GeoNormalizeVector (edgeHat);

			      if (edgeHat.y > VSMALL)
				{
				  visitedPoint[faces[face].node[j]] = 1;
				  nVisited++;

				  streamFunction[faces[face].node[j]] =
				    currentBStream + V_GetCmp (&uf,
							       face +
							       1) *
				    faces[face].Aj * LSGN (faces[face].n.x);
				}
			      else if (edgeHat.y < -VSMALL)
				{
				  visitedPoint[faces[face].node[j]] = 1;
				  nVisited++;

				  streamFunction[faces[face].node[j]] =
				    currentBStream - V_GetCmp (&uf,
							       face +
							       1) *
				    faces[face].Aj * LSGN (faces[face].n.x);
				}
			      else
				{
				  if (edgeHat.x > VSMALL)
				    {
				      visitedPoint[faces[face].node[j]] = 1;
				      nVisited++;

				      streamFunction[faces[face].node[j]] =
					currentBStream + V_GetCmp (&uf,
								   face +
								   1) *
					faces[face].Aj *
					LSGN (faces[face].n.y);
				    }
				  else if (edgeHat.x < -VSMALL)
				    {
				      visitedPoint[faces[face].node[j]] = 1;
				      nVisited++;

				      streamFunction[faces[face].node[j]] =
					currentBStream - V_GetCmp (&uf,
								   face +
								   1) *
					faces[face].Aj *
					LSGN (faces[face].n.y);
				    }
				}
			    }
			}

		    }


		}
	      else
		{
		  finished = LOGICAL_FALSE;
		}

	    }

	  // Internal faces

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

	      face = i;

	      pair = faces[face].pair;

	      if (pair == -1)
		continue;

	      pointFound = LOGICAL_FALSE;

	      currentStream = 0.0;

	      currentStreamPoint.x = 0.0;
	      currentStreamPoint.y = 0.0;
	      currentStreamPoint.z = 0.0;

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

		  // Check if the point has been visited
		  if (visitedPoint[faces[face].node[j]] == 1)
		    {
		      // The point has been visited
		      currentStream = streamFunction[faces[face].node[j]];
		      currentStreamPoint = nodes[faces[face].node[j]];

		      pointFound = LOGICAL_TRUE;

		      break;
		    }

		}

	      if (pointFound == LOGICAL_TRUE)
		{

		  // Sort out other points on the face                            
		  for (j = 0; j < faces[face].nbnodes; j++)
		    {

		      // Check if the point has been visited
		      if (visitedPoint[faces[face].node[j]] == 0)
			{

			  edgeHat =
			    GeoSubVectorVector (nodes[faces[face].node[j]],
						currentStreamPoint);

			  edgeHat.z = 0.0;

			  edgeHat = GeoNormalizeVector (edgeHat);

			  if (edgeHat.y > VSMALL)
			    {
			      visitedPoint[faces[face].node[j]] = 1;
			      nVisited++;

			      streamFunction[faces[face].node[j]] =
				currentStream + V_GetCmp (&uf,
							  face +
							  1) *
				faces[face].Aj * LSGN (faces[face].n.x);
			    }
			  else if (edgeHat.y < -VSMALL)
			    {
			      visitedPoint[faces[face].node[j]] = 1;
			      nVisited++;

			      streamFunction[faces[face].node[j]] =
				currentStream - V_GetCmp (&uf,
							  face +
							  1) *
				faces[face].Aj * LSGN (faces[face].n.x);
			    }
			}

		    }


		}
	      else
		{
		  finished = LOGICAL_FALSE;
		}


	    }

	  if (nVisited == nVisitedOld)
	    {

	      //printf("Exhausted a seed. Looking for new seed.\n");

	      break;
	    }
	  else
	    {
	      nVisitedOld = nVisited;
	    }

	}
      while (finished == LOGICAL_FALSE);


    }
  while (finished == LOGICAL_FALSE);

  free (visitedPoint);

  // Get maximum and minimum values       
  vmin = +VGREAT;
  vmax = -VGREAT;

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

      vmin = LMIN (vmin, streamFunction[i]);
      vmax = LMAX (vmax, streamFunction[i]);

    }

  // Normalize stream function

  if (vmax != 0.0)
    {
      for (i = 0; i < nbnodes; i++)
	{

	  streamFunction[i] = (streamFunction[i] - vmin) / (vmax - vmin);
	  //streamFunction[i] /= vmax; 

	}
    }

}

void
WriteProbeViews (FILE * fp, char *var, double curtime)
{

  int i, j, k;

  int node, element;

  double vs;
  double sd;

  double *se, *sf;

  vs = 0.0;
  
  se = calloc (nbelements, sizeof (double));
  sf = calloc (nbfaces, sizeof (double));

  // Create probe views
  for (k = 0; k < nphi; k++)
    {

      if (parameter.probe[k] == LOGICAL_FALSE)
	continue;

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

	  element = i;

	  if (k == 0)
	    vs = V_GetCmp (&xu, element + 1);
	  if (k == 1)
	    vs = V_GetCmp (&xv, element + 1);
	  if (k == 2)
	    vs = V_GetCmp (&xw, element + 1);
	  if (k == 3)
	    vs = V_GetCmp (&xp, element + 1);
	  if (k == 4)
	    vs = V_GetCmp (&xT, element + 1);
	  if (k == 5)
	    vs = V_GetCmp (&xs, element + 1);

	  se[element] = vs;

	}

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

         face = i;

         if (k == 0)
         vs = V_GetCmp (&xuf, face + 1);
         if (k == 1)
         vs = V_GetCmp (&xvf, face + 1);
         if (k == 2)
         vs = V_GetCmp (&xwf, face + 1);
         if (k == 3)
         vs = V_GetCmp (&xpf, face + 1);
         if (k == 4)
         vs = V_GetCmp (&xTf, face + 1);
         if (k == 5)
         vs = V_GetCmp (&xsf, face + 1);

         sf[face] = vs;

         }
       */

      fprintf (fp, "View \"Variable: %c\"{\n", var[k]);

      fprintf (fp, "TIME { %f };\n", curtime);

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

         face = i;

         if (faces[face].type == TRIANGLE) fprintf (fp, "ST(");
         if (faces[face].type == QUADRANGLE) fprintf (fp, "SQ(");

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

         fprintf (fp, "%f, %f, %f", nodes[node].x, 
         nodes[node].y, 
         nodes[node].z;

         if (j != faces[face].nbnodes - 1)
         fprintf (fp, ",");

         }

         fprintf (fp, ")");

         fprintf (fp, "{");

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

         sd = sf[face];

         fprintf (fp, " %f", sd);

         if (j != faces[face].nbnodes - 1)
         fprintf (fp, ",");

         }

         fprintf (fp, "};\n");

         }
       */

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

	  element = i;

	  if (elements[element].type == TETRAHEDRON)
	    fprintf (fp, "SS(");
	  if (elements[element].type == HEXAHEDRON)
	    fprintf (fp, "SH(");
	  if (elements[element].type == PRISM)
	    fprintf (fp, "SI(");

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

	      fprintf (fp, "%f, %f, %f", nodes[node].x,
		       nodes[node].y, nodes[node].z);

	      if (j != elements[element].nbnodes - 1)
		fprintf (fp, ",");

	    }