pos_print.c

cfd求解器使用与gmsh网格的求解

#define RCSID "$Id: Pos_Print.c,v 1.81 2006/02/26 00:42:59 geuzaine Exp $"
/*
 * Copyright (C) 1997-2006 P. Dular, C. Geuzaine
 *
 * 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.
 *
 * Please report all bugs and problems to <[e-mail]>.
 */

#include "GetDP.h"
#include "Data_Passive.h"
#include "Numeric.h"
#include "Treatment_Formulation.h"
#include "Get_Geometry.h"
#include "CurrentData.h"
#include "Get_DofOfElement.h"
#include "ExtendedGroup.h"
#include "Cal_Quantity.h"
#include "Pos_Formulation.h"
#include "Pos_Quantity.h"
#include "Pos_Element.h"
#include "Pos_Search.h"
#include "Pos_Print.h"
#include "Pos_Format.h"
#include "Adapt.h"
#include "Tools.h"

extern int  InteractiveInterrupt ;

/*
  Print OnElementsOf
  ------------------
  expl: plot on elements, belonging to the current mesh, where 
        the solution was computed during the processing stage
  args: list of groups of region type
  
  Print OnSection 
  ---------------
  expl: plot an a 'real' cut of the mesh, i.e. computation on the 
        intersections of the mesh with a cutting entity (plane, line)
  args: 2 (not done) or 3 points, specifying the cutting line or the cutting plane
  
  Print OnGrid
  ------------
  expl: reinterpolate the solution on a grid
  args: - a list of groups of region type (belonging to a mesh, where the
	  solution will be reinterpolated)
        - 3 expressions (using $S and $T) and 2 intervals for the parametric 
	  grid definition

  Print OnPoint, OnLine, OnPlane, OnBox
  -------------------------------------
  expl: reinterpolate the solution on a grid (particular cases)
  args: 1, 2, 3 or 4 points (0d, 1d, 2d or 3d grid) and the associated
        number of divisions

  Print OnRegion
  --------------
  expl: print Global Quantities associated with Regions
  args: list of groups of region type

*/

/* ------------------------------------------------------------------------ */
/*  P o s _ P r i n t O n E l e m e n t s O f                               */
/* ------------------------------------------------------------------------ */

struct CutEdge {
  int     nbc ;
  double  x[2],y[2],z[2] ;  
  double  xc,yc,zc ;
  double  uc,vc,wc ;
  struct  Value  *Value ;
} ;

struct xyzv{
  double x,y,z;
  struct Value v;
  /*int nbvals; for time domain -> malloc Value *v... */
  int nboccurences;
};

static int fcmp_xyzv(const void * a, const void * b) {
  struct xyzv *p1, *p2;
  double TOL=Current.GeoData->CharacteristicLength * 1.e-8;
  p1 = (struct xyzv*)a;
  p2 = (struct xyzv*)b;
  if(p1->x - p2->x > TOL) return 1;
  if(p1->x - p2->x <-TOL) return -1;
  if(p1->y - p2->y > TOL) return 1;
  if(p1->y - p2->y <-TOL) return -1;
  if(p1->z - p2->z > TOL) return 1;
  if(p1->z - p2->z <-TOL) return -1;
  return 0;
}

static List_T * SkinPostElement_L ;
static int      SkinDepth ;

static void Cut_SkinPostElement(void *a, void *b){
  struct PostElement  * PE ;

  PE = *(struct PostElement**)a ;

  Cut_PostElement(PE, Geo_GetGeoElement(PE->Index), SkinPostElement_L, 
		  PE->Index, SkinDepth, 0, 1) ;
}

static void Decompose_SkinPostElement(void *a, void *b){
  struct PostElement  * PE, * PE2 ;

  PE = *(struct PostElement**)a ;

  if(PE->Type != QUADRANGLE) return;
  /* change the quad to a tri */
  PE->Type = TRIANGLE;
  PE->NbrNodes = 3;
  /* create a second tri */
  PE2 = NodeCopy_PostElement(PE) ;
  PE2->NumNodes[1] = PE->NumNodes[2];
  PE2->u[1] = PE->u[2]; PE2->x[1] = PE->x[2];
  PE2->v[1] = PE->v[2]; PE2->y[1] = PE->y[2];
  PE2->w[1] = PE->w[2]; PE2->z[1] = PE->z[2];
  PE2->NumNodes[2] = PE->NumNodes[3];
  PE2->u[2] = PE->u[3]; PE2->x[2] = PE->x[3];
  PE2->v[2] = PE->v[3]; PE2->y[2] = PE->y[3];
  PE2->w[2] = PE->w[3]; PE2->z[2] = PE->z[3];
  List_Add(SkinPostElement_L, &PE2);
}

void  Pos_PrintOnElementsOf(struct PostQuantity     *NCPQ_P,
			    struct PostQuantity     *CPQ_P,
			    int                      Order,
			    struct DefineQuantity   *DefineQuantity_P0,
			    struct QuantityStorage  *QuantityStorage_P0,
			    struct PostSubOperation *PostSubOperation_P) {

  Tree_T  * PostElement_T ;
  List_T  * PostElement_L, * Region_L ;

  struct Element        Element ;
  struct PostElement  * PE ;
  struct Value        * CumulativeValues ;
  struct xyzv           xyzv, *xyzv_P ;
  Tree_T              * xyzv_T ;
  double  * Error=NULL, Dummy[5], d, x1, x2 ;
  int       jj, NbrGeo, iGeo, incGeo, NbrPost=0, iPost ;
  int       NbrTimeStep, iTime, iNode ;
  int       Store = 0, DecomposeInSimplex = 0, Depth ;

  GetDP_Begin("Pos_PrintOnElementsOf");

  /* Select the TimeSteps */
  
  NbrTimeStep = Pos_InitTimeSteps(PostSubOperation_P);
  
  /* Print the header */

  NbrGeo = Geo_GetNbrGeoElements() ;

  Format_PostHeader(PostSubOperation_P->Format, 
		    PostSubOperation_P->Iso, NbrTimeStep,
		    PostSubOperation_P->HarmonicToTime, 
		    PostSubOperation_P->CombinationType, Order,
		    NCPQ_P?NCPQ_P->Name:NULL,
		    CPQ_P?CPQ_P->Name:NULL);
  
  /* Get the region */
  
  Region_L = ((struct Group *)
	      List_Pointer(Problem_S.Group, 
			   PostSubOperation_P->Case.OnRegion.RegionIndex))->InitialList ;
  Get_InitDofOfElement(&Element) ;

  /* Compute the Cumulative quantity, if any */

  if(CPQ_P){
    Cal_PostCumulativeQuantity(Region_L,
			       PostSubOperation_P->PostQuantitySupport[Order],
			       PostSubOperation_P->TimeStep_L,
			       CPQ_P, DefineQuantity_P0,
			       QuantityStorage_P0, &CumulativeValues);
  }

  /* Init a search grid if we plot a NonCumulative quantity with OnGrid */

  if(NCPQ_P && PostSubOperation_P->SubType == PRINT_ONGRID)
    Init_SearchGrid(&Current.GeoData->Grid) ;

  /* If we compute a skin, apply smoothing, sort the results, or
     perform adaption, we'll need to store all the PostElements */

  if(PostSubOperation_P->Smoothing || PostSubOperation_P->Skin ||
     PostSubOperation_P->Adapt || PostSubOperation_P->Sort)
    Store = 1 ;

  /* Check if everything is OK for Adaption */

  if(PostSubOperation_P->Adapt){
    if(PostSubOperation_P->Dimension == _ALL)
      Msg(GERROR, "You have to specify a Dimension for the adaptation (2 or 3)");
    if(PostSubOperation_P->Target < 0.)
      Msg(GERROR, "You have to specify a Target for the adaptation (e.g. 0.01)");
    if(NbrTimeStep > 1)
      Msg(GERROR, "Adaption not ready with more than one time step"); 
  }

  /* Check if we should decompose all PostElements to simplices */

  if(!PostSubOperation_P->Skin && PostSubOperation_P->DecomposeInSimplex)
    DecomposeInSimplex = 1 ;

  /* Check for de-refinement */
  
  if(PostSubOperation_P->Depth < 0)
    incGeo = - PostSubOperation_P->Depth ;
  else
    incGeo = 1 ;
  
  /* Create the list of PostElements */

  PostElement_L = List_Create(Store ? NbrGeo/10 : 10, Store ? NbrGeo/10 : 10,
			      sizeof(struct PostElement *)) ;

  if(Store){
    
    /* If we have a Skin, we will divide after the skin extraction */

    if(PostSubOperation_P->Skin && PostSubOperation_P->Depth > 1)
      Depth = 1;
    else
      Depth = PostSubOperation_P->Depth;
    
    /* Generate all PostElements */

    for(iGeo = 0 ; iGeo < NbrGeo ; iGeo += incGeo) {
      if(InteractiveInterrupt) break;
      Progress(iGeo, NbrGeo, "Generate: ") ;
      Element.GeoElement = Geo_GetGeoElement(iGeo) ;
      if(List_Search(Region_L, &Element.GeoElement->Region, fcmp_int)){
	Fill_PostElement(Element.GeoElement, PostElement_L, iGeo,
			 Depth, PostSubOperation_P->Skin,
			 PostSubOperation_P->EvaluationPoints,
			 DecomposeInSimplex) ;
      }
    }

    /* Compute the skin */

    if(PostSubOperation_P->Skin){
      PostElement_T = Tree_Create(sizeof(struct PostElement *), fcmp_PostElement);

      for(iPost = 0 ; iPost < List_Nbr(PostElement_L) ; iPost++){
	if(InteractiveInterrupt) break;
	Progress(iPost, List_Nbr(PostElement_L), "Skin: ") ;
	PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
	if(Tree_PQuery(PostElement_T, &PE)){
	  Tree_Suppress(PostElement_T, &PE);
	  Destroy_PostElement(PE) ;
	}
	else
	  Tree_Add(PostElement_T, &PE);
      }

      /* only decompose in simplices (triangles!) now */
      if(PostSubOperation_P->DecomposeInSimplex){
	List_Reset(PostElement_L);
	SkinPostElement_L = PostElement_L ;
	Tree_Action(PostElement_T, Decompose_SkinPostElement);
	for(iPost = 0 ; iPost < List_Nbr(SkinPostElement_L) ; iPost++)
	  Tree_Add(PostElement_T, (struct PostElement**)List_Pointer(SkinPostElement_L, iPost)) ;
      }

      if(PostSubOperation_P->Depth > 1){
	List_Reset(PostElement_L);
	SkinPostElement_L = PostElement_L ;
	SkinDepth = PostSubOperation_P->Depth ;
	Tree_Action(PostElement_T, Cut_SkinPostElement) ;
      }
      else{
	List_Delete(PostElement_L);
	PostElement_L = Tree2List(PostElement_T);
      }

      Tree_Delete(PostElement_T);
    }

  } /* if Store */
  
  /* Loop on GeoElements */
  
  for(iGeo = 0 ; iGeo < NbrGeo ; iGeo += incGeo) {
    
    if(InteractiveInterrupt) break;
    
    if(Store){
      if(iGeo) break ;
    }
    else{
      Progress(iGeo, NbrGeo, "Compute: ") ;
      List_Reset(PostElement_L) ;
      Element.GeoElement = Geo_GetGeoElement(iGeo) ;
      if(List_Search(Region_L, &Element.GeoElement->Region, fcmp_int)){
	Fill_PostElement(Element.GeoElement, PostElement_L, iGeo,
			 PostSubOperation_P->Depth, PostSubOperation_P->Skin,
			 PostSubOperation_P->EvaluationPoints,
			 DecomposeInSimplex) ;
      }
    }

    NbrPost = List_Nbr(PostElement_L) ;  

    /* Loop on PostElements */
    
    for(iPost = 0 ; iPost < NbrPost ; iPost++) {

      if(InteractiveInterrupt) break;
      
      if(Store)
	Progress(iPost, NbrPost, "Compute: ") ;

      PE = *(struct PostElement **)List_Pointer(PostElement_L, iPost) ;

      if(!NCPQ_P){ /* Only one Cumulative */
	for (iTime = 0 ; iTime < NbrTimeStep ; iTime++){
	  for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++)
	    Cal_CopyValue(&CumulativeValues[iTime], &PE->Value[iNode]);
	  if(!Store)
	    Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso, 0,
			       Current.Time, iTime, NbrTimeStep, 
			       Current.NbrHar, PostSubOperation_P->HarmonicToTime,
			       NULL, PE, 
			       PostSubOperation_P->ChangeOfCoordinates,
			       PostSubOperation_P->ChangeOfValues);
	}
      }
      else{ /* There is one non-cumulative */

	if(PostSubOperation_P->SubType == PRINT_ONGRID){ /* We re-interpolate */
	  for (iTime = 0 ; iTime < NbrTimeStep ; iTime++) {
	    Pos_InitAllSolutions(PostSubOperation_P->TimeStep_L, iTime) ;
	    for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
	      InWhichElement(Current.GeoData->Grid, Region_L, &Element, 
			     PostSubOperation_P->Dimension, 
			     PE->x[iNode], PE->y[iNode], PE->z[iNode],
			     &PE->u[iNode], &PE->v[iNode], &PE->w[iNode]) ;
	      Current.Region = Element.Region ;
	      Current.x = PE->x[iNode] ; 
	      Current.y = PE->y[iNode] ;
	      Current.z = PE->z[iNode] ;
	      Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,
			       NULL, &Element,
			       PE->u[iNode], PE->v[iNode], PE->w[iNode], &PE->Value[iNode]);
	      if(CPQ_P)
		Combine_PostQuantity(PostSubOperation_P->CombinationType, Order,
				     &PE->Value[iNode], &CumulativeValues[iNode]) ;
	    }
	    if(!Store)
	      Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso, 0,
				 Current.Time, iTime, NbrTimeStep, 
				 Current.NbrHar, PostSubOperation_P->HarmonicToTime,
				 NULL, PE, 
				 PostSubOperation_P->ChangeOfCoordinates,
				 PostSubOperation_P->ChangeOfValues);
	  }
	}
	else{ /* PRINT_ONREGION: We work on the real mesh */
	  Element.GeoElement = Geo_GetGeoElement(PE->Index) ;
	  Element.Num = Element.GeoElement->Num ;
	  Element.Type = Element.GeoElement->Type ;
	  Current.Region = Element.Region = Element.GeoElement->Region ;
	  Get_NodesCoordinatesOfElement(&Element) ;

	  for (iTime = 0 ; iTime < NbrTimeStep ; iTime++) {

	    Pos_InitAllSolutions(PostSubOperation_P->TimeStep_L, iTime) ;

	    for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
	      Current.x = PE->x[iNode] ;
	      Current.y = PE->y[iNode] ;