#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 <getdp@geuz.org>.
*/
#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] ;
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[iTime]) ;
}
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);
}
}
}
if(!Store) Destroy_PostElement(PE) ;
}
} /* for iGeo */
/* Perform Smoothing */
if(PostSubOperation_P->Smoothing && !InteractiveInterrupt){
Msg(INFO, "Smoothing (phase 1)");
xyzv_T = Tree_Create(sizeof(struct xyzv), fcmp_xyzv);
for (iPost = 0 ; iPost < NbrPost ; iPost++){
PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++) {
xyzv.x = PE->x[iNode];
xyzv.y = PE->y[iNode];
xyzv.z = PE->z[iNode];
if((xyzv_P = (struct xyzv*)Tree_PQuery(xyzv_T, &xyzv))){
x1 = (double)(xyzv_P->nboccurences)/ (double)(xyzv_P->nboccurences + 1.);
x2 = 1./(double)(xyzv_P->nboccurences + 1);
Cal_AddMultValue2(&xyzv_P->v, x1, &PE->Value[iNode], x2);
xyzv_P->nboccurences++;
}
else{
Cal_CopyValue(&PE->Value[iNode],&xyzv.v);
xyzv.nboccurences = 1;
Tree_Add(xyzv_T, &xyzv);
}
}
}
Msg(INFO, "Smoothing (phase 2)");
for(iPost = 0 ; iPost < NbrPost ; iPost++){
PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost) ;
for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++) {
xyzv.x = PE->x[iNode];
xyzv.y = PE->y[iNode];
xyzv.z = PE->z[iNode];
if((xyzv_P = (struct xyzv*)Tree_PQuery(xyzv_T, &xyzv))){
Cal_CopyValue(&xyzv_P->v, &PE->Value[iNode]);
}
else
Msg(WARNING, "Node (%g,%g,%g) not found", xyzv.x, xyzv.y, xyzv.z);
}
}
Tree_Delete(xyzv_T);
} /* if Smoothing */
/* Perform Adaption */
if(PostSubOperation_P->Adapt && !InteractiveInterrupt){
if(!Current.GeoData->H)
Current.GeoData->H = (double*)Malloc((NbrGeo+2)*sizeof(double)) ;
if(!Current.GeoData->P)
Current.GeoData->P = (double*)Malloc((NbrGeo+2)*sizeof(double)) ;
Error = (double*)Malloc((NbrGeo+1)*sizeof(double)) ;
/* All elements are perfect... */
for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++){
Element.GeoElement = Geo_GetGeoElement(iGeo) ;
Element.Num = Element.GeoElement->Num ;
Element.Type = Element.GeoElement->Type ;
Element.Region = Element.GeoElement->Region ;
Get_NodesCoordinatesOfElement(&Element) ;
Current.GeoData->H[iGeo+1] = Cal_MaxEdgeLength(&Element) ;
Current.GeoData->P[iGeo+1] = 1. ;
Error[iGeo+1] = PostSubOperation_P->Target ;
}
/* ...except those we want to optimize */
for(iPost = 0 ; iPost < NbrPost ; iPost++){
PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
Error[PE->Index+1] = 0. ;
for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++)
Error[PE->Index+1] += PE->Value[iNode].Val[0] ;
Error[PE->Index+1] /= (double)PE->NbrNodes ;
}
Adapt (NbrGeo, PostSubOperation_P->Adapt,
PostSubOperation_P->Dimension, Error,
Current.GeoData->H, Current.GeoData->P,
PostSubOperation_P->Target);
/* Clean up the interpolation orders to fit to what's available */
if(List_Nbr(PostSubOperation_P->Value_L)){
for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++){
for(jj = List_Nbr(PostSubOperation_P->Value_L)-1 ; jj >= 0 ; jj--){
d = *(double*)List_Pointer(PostSubOperation_P->Value_L, jj);
if(Current.GeoData->P[iGeo+1] > d || jj == 0){
Current.GeoData->P[iGeo+1] = d ;
break ;
}
}
}
}
} /* if Adapt */
/* Print everything if we are in Store mode */
if(Store && !InteractiveInterrupt){
/* Sort the elements */
switch(PostSubOperation_P->Sort){
case SORT_BY_POSITION : List_Sort(PostElement_L, fcmp_PostElement) ; break ;
case SORT_BY_CONNECTIVITY : Sort_PostElement_Connectivity(PostElement_L) ; break ;
}
Dummy[0] = Dummy[1] = Dummy[2] = Dummy[3] = Dummy[4] = 0. ;
for(iPost = 0 ; iPost < NbrPost ; iPost++){
PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
/* Get the values from adaption */
if(PostSubOperation_P->Adapt){
Element.GeoElement = Geo_GetGeoElement(PE->Index) ;
Dummy[0] = Element.GeoElement->Num ;
Dummy[1] = Error[PE->Index+1] ;
Dummy[2] = Current.GeoData->H[PE->Index+1] ;
Dummy[3] = Current.GeoData->P[PE->Index+1] ;
Dummy[4] = iPost ? 0 : NbrPost ;
for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){
PE->Value[iNode].Type = SCALAR ;
if(PostSubOperation_P->Adapt == H1 ||
PostSubOperation_P->Adapt == H2)
PE->Value[iNode].Val[0] = Dummy[2] ;
else
PE->Value[iNode].Val[0] = Dummy[3] ;
}
}
/* Compute curvilinear coord if connection sort */
if(PostSubOperation_P->Sort == SORT_BY_CONNECTIVITY){
Dummy[0] = Dummy[1] ;
Dummy[1] = Dummy[0] + sqrt(DSQU(PE->x[1]-PE->x[0])+
DSQU(PE->y[1]-PE->y[0])+
DSQU(PE->z[1]-PE->z[0])) ;
Dummy[2] = PE->v[0] ;
Dummy[3] = -1. ;
}
Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso, 1,
Current.Time, 0, 1,
Current.NbrHar, PostSubOperation_P->HarmonicToTime,
Dummy, PE,
PostSubOperation_P->ChangeOfCoordinates,
PostSubOperation_P->ChangeOfValues);
}
}
Format_PostFooter(PostSubOperation_P, Store);
if(Store)
for(iPost = 0 ; iPost < NbrPost ; iPost++){
PE = *(struct PostElement**)List_Pointer(PostElement_L, iPost);
Destroy_PostElement(PE) ;
}
List_Delete(PostElement_L);
if(CPQ_P) Free(CumulativeValues);
if(PostSubOperation_P->Adapt) Free(Error) ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* P o s _ P r i n t O n S e c t i o n */
/* ------------------------------------------------------------------------ */
double Plane(double a, double b, double c, double d,
double x, double y, double z){
GetDP_Begin("Plane");
GetDP_Return(a*x+b*y+c*z+d);
}
static double DIRX[3], DIRY[3], DIRZ[3], XCP, YCP ;
int fcmp_Angle (const void *a, const void *b){
struct CutEdge *q , *w;
double x1,y1,x2,y2,ang1,ang2;
q = (struct CutEdge*)a;
w = (struct CutEdge*)b;
x1 = q->xc*DIRX[0] + q->yc*DIRX[1] + q->zc*DIRX[2];
y1 = q->xc*DIRY[0] + q->yc*DIRY[1] + q->zc*DIRY[2];
x2 = w->xc*DIRX[0] + w->yc*DIRX[1] + w->zc*DIRX[2];
y2 = w->xc*DIRY[0] + w->yc*DIRY[1] + w->zc*DIRY[2];
ang1 = atan2(y1-YCP,x1-XCP);
ang2 = atan2(y2-YCP,x2-XCP);
if(ang1>ang2)return 1;
return -1;
}
void prodvec (double *a , double *b , double *c){
GetDP_Begin("prodvec");
c[0] = a[1]*b[2]-a[2]*b[1];
c[1] = a[2]*b[0]-a[0]*b[2];
c[2] = a[0]*b[1]-a[1]*b[0];
GetDP_End ;
}
void normvec(double *a){
double mod;
GetDP_Begin("normvec");
mod = sqrt(SQU(a[0])+SQU(a[1])+SQU(a[2]));
a[0]/=mod;
a[1]/=mod;
a[2]/=mod;
GetDP_End ;
}
#define NBR_MAX_CUT 10
#define LETS_PRINT_THE_RESULT \
List_Reset(PE_L); \
if(PostSubOperation_P->Depth < 2) \
List_Add(PE_L, &PE) ; \
else \
Cut_PostElement(PE, Element.GeoElement, PE_L, PE->Index, \
PostSubOperation_P->Depth, 0, 1) ; \
for(iPost = 0 ; iPost < List_Nbr(PE_L) ; iPost++){ \
PE = *(struct PostElement **)List_Pointer(PE_L, iPost) ; \
for(iTime = 0 ; iTime < NbTimeStep ; iTime++){ \
Pos_InitAllSolutions(PostSubOperation_P->TimeStep_L, iTime) ; \
for(iNode = 0 ; iNode < PE->NbrNodes ; iNode++){ \
if(NCPQ_P){ \
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[iTime]) ;\
} \
else \
Cal_CopyValue(&CumulativeValues[iTime],&PE->Value[iNode]); \
} \
Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso,0, \
Current.Time, iTime, NbTimeStep, \
Current.NbrHar, PostSubOperation_P->HarmonicToTime, \
NULL, PE, \
PostSubOperation_P->ChangeOfCoordinates, \
PostSubOperation_P->ChangeOfValues); \
} \
} \
for(iPost = 0 ; iPost < List_Nbr(PE_L) ; iPost++) \
Destroy_PostElement(*(struct PostElement **)List_Pointer(PE_L, iPost));
void Pos_PrintOnSection(struct PostQuantity *NCPQ_P,
struct PostQuantity *CPQ_P,
int Order,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct PostSubOperation *PostSubOperation_P) {
struct CutEdge e[NBR_MAX_CUT];
struct Element Element ;
struct PostElement * PE ;
struct Value * CumulativeValues ;
List_T * PE_L ;
int NbGeoElement, NbTimeStep, NbCut, * NumNodes ;
int iPost, iNode, iGeo, iCut, iEdge, iTime ;
double A, B, C, D, d1, d2, u, xcg, ycg, zcg ;
double x[3], y[3], z[3] ;
GetDP_Begin("Pos_PrintOnSection");
NbTimeStep = Pos_InitTimeSteps(PostSubOperation_P);
PE_L = List_Create(10, 10, sizeof(struct PostElement *)) ;
for(iCut = 0 ; iCut < NBR_MAX_CUT ; iCut++)
e[iCut].Value = (struct Value*) Malloc(NbTimeStep*sizeof(struct Value)) ;
Format_PostHeader(PostSubOperation_P->Format,
PostSubOperation_P->Iso, NbTimeStep,
PostSubOperation_P->HarmonicToTime,
PostSubOperation_P->CombinationType, Order,
NCPQ_P?NCPQ_P->Name:NULL,
CPQ_P?CPQ_P->Name:NULL);
if(CPQ_P){
Cal_PostCumulativeQuantity(NULL,
PostSubOperation_P->PostQuantitySupport[Order],
PostSubOperation_P->TimeStep_L,
CPQ_P, DefineQuantity_P0,
QuantityStorage_P0, &CumulativeValues);
}
switch(PostSubOperation_P->SubType) {
case PRINT_ONSECTION_1D :
Msg(GERROR, "Print on 1D cuts not done (use Print OnLine instead)");
break;
case PRINT_ONSECTION_2D :
/* Ax+By+Cz+D=0 from (x0,y0,z0),(x1,y1,z1),(x2,y2,z2) */
x[0] = PostSubOperation_P->Case.OnSection.x[0] ;
y[0] = PostSubOperation_P->Case.OnSection.y[0] ;
z[0] = PostSubOperation_P->Case.OnSection.z[0] ;
x[1] = PostSubOperation_P->Case.OnSection.x[1] ;
y[1] = PostSubOperation_P->Case.OnSection.y[1] ;
z[1] = PostSubOperation_P->Case.OnSection.z[1] ;
x[2] = PostSubOperation_P->Case.OnSection.x[2] ;
y[2] = PostSubOperation_P->Case.OnSection.y[2] ;
z[2] = PostSubOperation_P->Case.OnSection.z[2] ;
A = (y[1]-y[0])*(z[2]-z[0]) - (z[1]-z[0])*(y[2]-y[0]) ;
B = -(x[1]-x[0])*(z[2]-z[0]) + (z[1]-z[0])*(x[2]-x[0]) ;
C = (x[1]-x[0])*(y[2]-y[0]) - (y[1]-y[0])*(x[2]-x[0]) ;
D = -A*x[0]-B*y[0]-C*z[0] ;
/* Cut each element */
NbGeoElement = Geo_GetNbrGeoElements() ;
for(iGeo = 0 ; iGeo < NbGeoElement ; iGeo++) {
if(InteractiveInterrupt) break;
Progress(iGeo, NbGeoElement, "Cut: ") ;
Element.GeoElement = Geo_GetGeoElement(iGeo) ;
Element.Num = Element.GeoElement->Num ;
Element.Type = Element.GeoElement->Type ;
Current.Region = Element.Region = Element.GeoElement->Region ;
if((PostSubOperation_P->Dimension == _ALL &&
(Element.GeoElement->Type != POINT)) ||
(PostSubOperation_P->Dimension == _3D &&
(Element.GeoElement->Type & (TETRAHEDRON|HEXAHEDRON|PRISM|PYRAMID))) ||
(PostSubOperation_P->Dimension == _2D &&
(Element.GeoElement->Type & (TRIANGLE|QUADRANGLE))) ||
(PostSubOperation_P->Dimension == _1D &&
(Element.GeoElement->Type & LINE))){
Get_NodesCoordinatesOfElement(&Element) ;
if(Element.GeoElement->NbrEdges == 0)
Geo_CreateEdgesOfElement(Element.GeoElement) ;
NbCut = 0;
for(iEdge = 0 ; iEdge < Element.GeoElement->NbrEdges ; iEdge++){
NumNodes = Geo_GetNodesOfEdgeInElement(Element.GeoElement, iEdge) ;
e[NbCut].x[0] = Element.x[abs(NumNodes[0])-1] ;
e[NbCut].y[0] = Element.y[abs(NumNodes[0])-1] ;
e[NbCut].z[0] = Element.z[abs(NumNodes[0])-1] ;
e[NbCut].x[1] = Element.x[abs(NumNodes[1])-1] ;
e[NbCut].y[1] = Element.y[abs(NumNodes[1])-1] ;
e[NbCut].z[1] = Element.z[abs(NumNodes[1])-1] ;
d1 = Plane(A,B,C,D,e[NbCut].x[0],e[NbCut].y[0],e[NbCut].z[0]);
d2 = Plane(A,B,C,D,e[NbCut].x[1],e[NbCut].y[1],e[NbCut].z[1]);
if(d1*d2 <= 0) {
if(d1*d2 < 0.) u = -d2/(d1-d2) ;
else if(d1 == 0.) u = 1. ;
else u = 0. ;
e[NbCut].xc = u*e[NbCut].x[0] + (1.-u)*e[NbCut].x[1];
e[NbCut].yc = u*e[NbCut].y[0] + (1.-u)*e[NbCut].y[1];
e[NbCut].zc = u*e[NbCut].z[0] + (1.-u)*e[NbCut].z[1];
if(NCPQ_P)
xyz2uvwInAnElement(&Element, e[NbCut].xc, e[NbCut].yc, e[NbCut].zc,
&e[NbCut].uc, &e[NbCut].vc, &e[NbCut].wc);
NbCut++;
}
}
if(NbCut > 3){
xcg = ycg = zcg = 0.;
for(iCut = 0 ; iCut < NbCut ; iCut++){
xcg += e[iCut].xc; ycg += e[iCut].yc; zcg += e[iCut].zc;
}
xcg /= (double)NbCut; ycg /= (double)NbCut; zcg /= (double)NbCut;
DIRZ[0] = A; DIRY[0] = xcg-e[0].xc;
DIRZ[1] = B; DIRY[1] = ycg-e[0].yc;
DIRZ[2] = C; DIRY[2] = zcg-e[0].zc;
normvec(DIRZ); normvec(DIRY); prodvec(DIRY,DIRZ,DIRX); normvec(DIRX);
XCP = xcg*DIRX[0] + ycg*DIRX[1] + zcg*DIRX[2];
YCP = xcg*DIRY[0] + ycg*DIRY[1] + zcg*DIRY[2];
qsort(e,NbCut,sizeof(struct CutEdge), fcmp_Angle);
}
if(NbCut > 2){
iCut = 2;
while(iCut < NbCut){
if(PostSubOperation_P->Depth > 0){
PE = Create_PostElement(iGeo, TRIANGLE, 3, 1) ;
PE->x[0] = e[0].xc; PE->x[1] = e[iCut-1].xc; PE->x[2] = e[iCut].xc;
PE->y[0] = e[0].yc; PE->y[1] = e[iCut-1].yc; PE->y[2] = e[iCut].yc;
PE->z[0] = e[0].zc; PE->z[1] = e[iCut-1].zc; PE->z[2] = e[iCut].zc;
PE->u[0] = e[0].uc; PE->u[1] = e[iCut-1].uc; PE->u[2] = e[iCut].uc;
PE->v[0] = e[0].vc; PE->v[1] = e[iCut-1].vc; PE->v[2] = e[iCut].vc;
PE->w[0] = e[0].wc; PE->w[1] = e[iCut-1].wc; PE->w[2] = e[iCut].wc;
LETS_PRINT_THE_RESULT ;
}
else{
PE = Create_PostElement(iGeo, POINT, 1, 0) ;
PE->x[0] = (e[0].xc + e[iCut-1].xc + e[iCut].xc) / 3. ;
PE->y[0] = (e[0].yc + e[iCut-1].yc + e[iCut].yc) / 3. ;
PE->z[0] = (e[0].zc + e[iCut-1].zc + e[iCut].zc) / 3. ;
PE->u[0] = (e[0].uc + e[iCut-1].uc + e[iCut].uc) / 3. ;
PE->v[0] = (e[0].vc + e[iCut-1].vc + e[iCut].vc) / 3. ;
PE->w[0] = (e[0].wc + e[iCut-1].wc + e[iCut].wc) / 3. ;
LETS_PRINT_THE_RESULT ;
}
iCut++;
}
}
if(NbCut == 2){
if(PostSubOperation_P->Depth > 0){
PE = Create_PostElement(iGeo, LINE, 2, 1) ;
PE->x[0] = e[0].xc; PE->x[1] = e[1].xc;
PE->y[0] = e[0].yc; PE->y[1] = e[1].yc;
PE->z[0] = e[0].zc; PE->z[1] = e[1].zc;
PE->u[0] = e[0].uc; PE->u[1] = e[1].uc;
PE->v[0] = e[0].vc; PE->v[1] = e[1].vc;
PE->w[0] = e[0].wc; PE->w[1] = e[1].wc;
LETS_PRINT_THE_RESULT ;
}
else{
PE = Create_PostElement(iGeo, POINT, 1, 0) ;
PE->x[0] = (e[0].xc + e[1].xc) / 2. ;
PE->y[0] = (e[0].yc + e[1].yc) / 2. ;
PE->z[0] = (e[0].zc + e[1].zc) / 2. ;
PE->u[0] = (e[0].uc + e[1].uc) / 2. ;
PE->v[0] = (e[0].vc + e[1].vc) / 2. ;
PE->w[0] = (e[0].wc + e[1].wc) / 2. ;
LETS_PRINT_THE_RESULT ;
}
}
}
}
Format_PostFooter(PostSubOperation_P, 0);
break;
default :
Msg(GERROR, "Unknown operation in Print OnSection");
break;
}
List_Delete(PE_L) ;
if(CPQ_P) Free(CumulativeValues);
for(iCut = 0 ; iCut < NBR_MAX_CUT ; iCut++) Free(e[iCut].Value) ;
GetDP_End ;
}
#undef NBR_MAX_CUT
#undef LETS_PRINT_THE_RESULT
/* ------------------------------------------------------------------------ */
/* P o s _ P r i n t O n G r i d */
/* ------------------------------------------------------------------------ */
#define LETS_PRINT_THE_RESULT \
PE->x[0] = Current.xp = Current.x ; \
PE->y[0] = Current.yp = Current.y ; \
PE->z[0] = Current.zp = Current.z ; \
if(!NCPQ_P){ \
for (ts = 0 ; ts < NbTimeStep ; ts++){ \
PE->Value[0] = CumulativeValues[ts] ; \
Format_PostElement(PSO_P->Format, PSO_P->Iso, 0, \
Current.Time, ts, NbTimeStep, \
Current.NbrHar, PSO_P->HarmonicToTime, \
Normal, PE, \
PSO_P->ChangeOfCoordinates, \
PSO_P->ChangeOfValues); \
} \
} \
else{ \
InWhichElement(Current.GeoData->Grid, NULL, &Element, PSO_P->Dimension, \
Current.x, Current.y, Current.z, &u, &v, &w) ; \
Current.Region = Element.Region ; \
if(Element.Num != NO_ELEMENT) \
PE->Index = Geo_GetGeoElementIndex(Element.GeoElement) ; \
else \
PE->Index = NO_ELEMENT ; \
for (ts = 0 ; ts < NbTimeStep ; ts++) { \
Pos_InitAllSolutions(PSO_P->TimeStep_L, ts) ; \
Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0, \
NULL, &Element, u, v, w, &PE->Value[0]); \
if(CPQ_P) \
Combine_PostQuantity(PSO_P->CombinationType, Order, \
&PE->Value[0], &CumulativeValues[ts]) ; \
Format_PostElement(PSO_P->Format, PSO_P->Iso, 0, \
Current.Time, ts, NbTimeStep, \
Current.NbrHar, PSO_P->HarmonicToTime, \
Normal, PE, \
PSO_P->ChangeOfCoordinates, \
PSO_P->ChangeOfValues); \
} \
}
#define ARRAY(i,j,k,t) \
Array[ (t) * Current.NbrHar * ((int)N[0]+1) * ((int)N[1]+1) + \
(k) * ((int)N[0]+1) * ((int)N[1]+1) + \
(j) * ((int)N[0]+1) + \
(i) ]
#define LETS_STORE_THE_RESULT \
if(!NCPQ_P){ \
if(CumulativeValues[0].Type != SCALAR) \
Msg(GERROR, "Print OnPlane is not designed for non scalar values with Depth > 1"); \
else \
for (ts = 0 ; ts < NbTimeStep ; ts++) \
for(k = 0 ; k < Current.NbrHar ; k++) \
ARRAY(i1,i2,k,ts) = (float)CumulativeValues[ts].Val[MAX_DIM*k] ; \
} \
else{ \
InWhichElement(Current.GeoData->Grid, NULL, &Element, PSO_P->Dimension, \
Current.x, Current.y, Current.z, &u, &v, &w) ; \
Current.Region = Element.Region ; \
for (ts = 0 ; ts < NbTimeStep ; ts++) { \
Pos_InitAllSolutions(PSO_P->TimeStep_L, ts) ; \
Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0, \
NULL, &Element, u, v, w, &PE->Value[0]); \
if(PE->Value[0].Type != SCALAR) \
Msg(GERROR, "Print OnPlane is not designed for non scalar values with Depth > 1");\
if(CPQ_P) \
Combine_PostQuantity(PSO_P->CombinationType, Order, \
&PE->Value[0], &CumulativeValues[ts]) ; \
for(k = 0 ; k < Current.NbrHar ; k++) \
ARRAY(i1,i2,k,ts) = (float)PE->Value[0].Val[MAX_DIM*k] ; \
} \
}
void Pos_PrintOnGrid(struct PostQuantity *NCPQ_P,
struct PostQuantity *CPQ_P,
int Order,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct PostSubOperation *PSO_P) {
struct Element Element ;
struct Value * CumulativeValues, Value ;
struct PostElement * PE , * PE2 ;
int i1, i2, i3, k, NbTimeStep, ts ;
float *Array = NULL ;
double u, v, w, Length, Normal[4] = {0., 0., 0., 0.} ;
double X[4], Y[4], Z[4], S[4], N[4], tmp[3];
GetDP_Begin("Pos_PrintOnGrid");
Get_InitDofOfElement(&Element) ;
NbTimeStep = Pos_InitTimeSteps(PSO_P);
if(CPQ_P){
Cal_PostCumulativeQuantity(NULL, PSO_P->PostQuantitySupport[Order],
PSO_P->TimeStep_L,
CPQ_P, DefineQuantity_P0,
QuantityStorage_P0, &CumulativeValues);
}
Init_SearchGrid(&Current.GeoData->Grid) ;
Format_PostHeader(PSO_P->Format, PSO_P->Iso,
NbTimeStep, PSO_P->HarmonicToTime,
PSO_P->CombinationType, Order,
NCPQ_P?NCPQ_P->Name:NULL,
CPQ_P?CPQ_P->Name:NULL);
PE = Create_PostElement(0, POINT, 1, 0) ;
switch(PSO_P->SubType) {
case PRINT_ONGRID_0D :
Current.x = PSO_P->Case.OnGrid.x[0] ;
Current.y = PSO_P->Case.OnGrid.y[0] ;
Current.z = PSO_P->Case.OnGrid.z[0] ;
Normal[0] = Normal[1] = Normal[2] = 0.0 ;
LETS_PRINT_THE_RESULT ;
break;
case PRINT_ONGRID_1D :
X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
N[0] = PSO_P->Case.OnGrid.n[0] ;
Normal[1] = Normal[2] = 0.0 ;
Length = sqrt(DSQU(X[1]-X[0]) + DSQU(Y[1]-Y[0]) + DSQU(Z[1]-Z[0])) ;
for (i1 = 0 ; i1 <= N[0] ; i1++) {
S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
Normal[0] = Length * S[0] ;
Current.x = X[0] + (X[1] - X[0]) * S[0] ;
Current.y = Y[0] + (Y[1] - Y[0]) * S[0] ;
Current.z = Z[0] + (Z[1] - Z[0]) * S[0] ;
LETS_PRINT_THE_RESULT ;
}
break;
case PRINT_ONGRID_2D :
X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
X[2] = PSO_P->Case.OnGrid.x[2] ;
Y[2] = PSO_P->Case.OnGrid.y[2] ;
Z[2] = PSO_P->Case.OnGrid.z[2] ;
S[0] = X[1]-X[0]; S[1] = Y[1]-Y[0]; S[2] = Z[1]-Z[0];
N[0] = X[2]-X[0]; N[1] = Y[2]-Y[0]; N[2] = Z[2]-Z[0];
prodvec(S,N,Normal);
Length = sqrt(DSQU(Normal[0])+DSQU(Normal[1])+DSQU(Normal[2]));
if(!Length){
Msg(WARNING, "Bad plane (null normal)");
GetDP_End ;
}
Normal[0]/=Length ; Normal[1]/=Length ; Normal[2]/=Length ;
N[0] = PSO_P->Case.OnGrid.n[0] ; N[1] = PSO_P->Case.OnGrid.n[1] ;
if(PSO_P->Depth > 1)
Array = (float*)
Malloc(NbTimeStep*Current.NbrHar*(int)((N[0]+1)*(N[1]+1))*sizeof(float)) ;
for (i1 = 0 ; i1 <= N[0] ; i1++) {
S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
for (i2 = 0 ; i2 <= N[1] ; i2++) {
S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
Current.x = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[1] ;
Current.y = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[1] ;
Current.z = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[1] ;
if(PSO_P->Depth > 1){
LETS_STORE_THE_RESULT ;
}
else{
LETS_PRINT_THE_RESULT ;
}
}
if(PSO_P->Depth < 2 && !Flag_BIN) fprintf(PostStream, "\n");
}
if(PSO_P->Depth > 1){
PE2 = Create_PostElement(0, TRIANGLE, 3, 0);
PE2->Value[0].Type = PE2->Value[1].Type = PE2->Value[2].Type = SCALAR ;
for (i1 = 0 ; i1 < N[0] ; i1++) {
S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
S[2] = (double)(i1+1) / (double)(N[0] ? N[0] : 1) ;
for (i2 = 0 ; i2 < N[1] ; i2++) {
S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
S[3] = (double)(i2+1) / (double)(N[1] ? N[1] : 1) ;
PE2->x[0] = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[1] ;
PE2->y[0] = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[1] ;
PE2->z[0] = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[1] ;
PE2->x[1] = X[0] + (X[1] - X[0]) * S[2] + (X[2] - X[0]) * S[1] ;
PE2->y[1] = Y[0] + (Y[1] - Y[0]) * S[2] + (Y[2] - Y[0]) * S[1] ;
PE2->z[1] = Z[0] + (Z[1] - Z[0]) * S[2] + (Z[2] - Z[0]) * S[1] ;
PE2->x[2] = X[0] + (X[1] - X[0]) * S[0] + (X[2] - X[0]) * S[3] ;
PE2->y[2] = Y[0] + (Y[1] - Y[0]) * S[0] + (Y[2] - Y[0]) * S[3] ;
PE2->z[2] = Z[0] + (Z[1] - Z[0]) * S[0] + (Z[2] - Z[0]) * S[3] ;
for (ts = 0 ; ts < NbTimeStep ; ts++){
for(k = 0 ; k < Current.NbrHar ; k++){
PE2->Value[0].Val[MAX_DIM*k] = ARRAY(i1,i2,k,ts) ;
PE2->Value[1].Val[MAX_DIM*k] = ARRAY(i1+1,i2,k,ts) ;
PE2->Value[2].Val[MAX_DIM*k] = ARRAY(i1,i2+1,k,ts) ;
}
Format_PostElement(PSO_P->Format, PSO_P->Iso, 0,
Current.Time, ts, NbTimeStep,
Current.NbrHar, PSO_P->HarmonicToTime,
Normal, PE2,
PSO_P->ChangeOfCoordinates,
PSO_P->ChangeOfValues);
}
PE2->x[0] = X[0] + (X[1] - X[0]) * S[2] + (X[2] - X[0]) * S[3] ;
PE2->y[0] = Y[0] + (Y[1] - Y[0]) * S[2] + (Y[2] - Y[0]) * S[3] ;
PE2->z[0] = Z[0] + (Z[1] - Z[0]) * S[2] + (Z[2] - Z[0]) * S[3] ;
tmp[0] = PE2->x[1]; PE2->x[1] = PE2->x[2]; PE2->x[2] = tmp[0];
tmp[1] = PE2->y[1]; PE2->y[1] = PE2->y[2]; PE2->y[2] = tmp[1];
tmp[2] = PE2->z[1]; PE2->z[1] = PE2->z[2]; PE2->z[2] = tmp[2];
for (ts = 0 ; ts < NbTimeStep ; ts++){
for(k = 0 ; k < Current.NbrHar ; k++){
PE2->Value[0].Val[MAX_DIM*k] = ARRAY(i1+1,i2+1,k,ts) ;
PE2->Value[1].Val[MAX_DIM*k] = ARRAY(i1,i2+1,k,ts) ;
PE2->Value[2].Val[MAX_DIM*k] = ARRAY(i1+1,i2,k,ts) ;
}
Format_PostElement(PSO_P->Format, PSO_P->Iso, 0,
Current.Time, ts, NbTimeStep,
Current.NbrHar, PSO_P->HarmonicToTime,
Normal, PE2,
PSO_P->ChangeOfCoordinates,
PSO_P->ChangeOfValues);
}
}
}
Destroy_PostElement(PE2) ;
Free(Array) ;
}
break;
case PRINT_ONGRID_3D :
X[0] = PSO_P->Case.OnGrid.x[0] ; X[1] = PSO_P->Case.OnGrid.x[1] ;
Y[0] = PSO_P->Case.OnGrid.y[0] ; Y[1] = PSO_P->Case.OnGrid.y[1] ;
Z[0] = PSO_P->Case.OnGrid.z[0] ; Z[1] = PSO_P->Case.OnGrid.z[1] ;
X[2] = PSO_P->Case.OnGrid.x[2] ; X[3] = PSO_P->Case.OnGrid.x[3] ;
Y[2] = PSO_P->Case.OnGrid.y[2] ; Y[3] = PSO_P->Case.OnGrid.y[3] ;
Z[2] = PSO_P->Case.OnGrid.z[2] ; Z[3] = PSO_P->Case.OnGrid.z[3] ;
N[0] = PSO_P->Case.OnGrid.n[0] ;
N[1] = PSO_P->Case.OnGrid.n[1] ;
N[2] = PSO_P->Case.OnGrid.n[2] ;
Normal[0] = Normal[1] = Normal[2] = 0.0 ;
for (i1 = 0 ; i1 <= N[0] ; i1++) {
S[0] = (double)i1 / (double)(N[0] ? N[0] : 1) ;
for (i2 = 0 ; i2 <= N[1] ; i2++) {
S[1] = (double)i2 / (double)(N[1] ? N[1] : 1) ;
for (i3 = 0 ; i3 <= N[2] ; i3++) {
S[2] = (double)i3 / (double)(N[2] ? N[2] : 1) ;
Current.x = X[0] + (X[1]-X[0])*S[0] + (X[2]-X[0])*S[1] + (X[3]-X[0])*S[2] ;
Current.y = Y[0] + (Y[1]-Y[0])*S[0] + (Y[2]-Y[0])*S[1] + (Y[3]-Y[0])*S[2] ;
Current.z = Z[0] + (Z[1]-Z[0])*S[0] + (Z[2]-Z[0])*S[1] + (Z[3]-Z[0])*S[2] ;
LETS_PRINT_THE_RESULT ;
}
if(!Flag_BIN) fprintf(PostStream, "\n");
}
if(!Flag_BIN) fprintf(PostStream, "\n\n");
/* two blanks lines for -index in gnuplot */
}
break;
case PRINT_ONGRID_PARAM :
for (i1 = 0 ; i1 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[0]) ; i1++) {
List_Read(PSO_P->Case.OnParamGrid.ParameterValue[0], i1, &Current.a) ;
for (i2 = 0 ; i2 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[1]) ; i2++) {
List_Read(PSO_P->Case.OnParamGrid.ParameterValue[1], i2, &Current.b) ;
for (i3 = 0 ; i3 < List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[2]) ; i3++) {
List_Read(PSO_P->Case.OnParamGrid.ParameterValue[2], i3, &Current.c) ;
Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[0],
NULL, 0., 0., 0., &Value) ;
Current.x = Value.Val[0];
Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[1],
NULL, 0., 0., 0., &Value) ;
Current.y = Value.Val[0];
Get_ValueOfExpressionByIndex(PSO_P->Case.OnParamGrid.ExpressionIndex[2],
NULL, 0., 0., 0., &Value) ;
Current.z = Value.Val[0];
Normal[0] = Current.a ;
Normal[1] = Current.b ;
Normal[2] = Current.c ;
LETS_PRINT_THE_RESULT ;
}
if(List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[2])>1 && !Flag_BIN)
fprintf(PostStream, "\n");
}
if(List_Nbr(PSO_P->Case.OnParamGrid.ParameterValue[1])>1 && !Flag_BIN)
fprintf(PostStream, "\n\n");
/* two blanks lines for -index in gnuplot */
}
break;
}
Destroy_PostElement(PE) ;
Format_PostFooter(PSO_P, 0);
if(CPQ_P) Free(CumulativeValues);
GetDP_End ;
}
#undef LETS_PRINT_THE_RESULT
#undef LETS_STORE_THE_RESULT
#undef ARRAY
/* ------------------------------------------------------------------------ */
/* P o s _ P r i n t O n R e g i o n */
/* ------------------------------------------------------------------------ */
void Pos_PrintOnRegion(struct PostQuantity *NCPQ_P,
struct PostQuantity *CPQ_P,
int Order,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct PostSubOperation *PostSubOperation_P) {
struct Element Element ;
struct Value Value, ValueSummed ;
struct PostQuantity *PQ_P ;
struct Group * Group_P ;
List_T *Region_L, *Support_L ;
int i, iTime, NbrTimeStep ;
int Nbr_Region=0, Num_Region, Group_FunctionType ;
int Flag_Summation=0, Type_Evaluation=0;
double u, v, w;
GetDP_Begin("Pos_PrintOnRegion");
NbrTimeStep = Pos_InitTimeSteps(PostSubOperation_P);
if (CPQ_P && NCPQ_P)
Msg(GERROR, "Only one PostProcessing Quantity allowed in PostOperation") ;
if (CPQ_P) {
PQ_P = CPQ_P ;
Support_L = /* for e.g. PQ[ Support ] ... */
((struct Group *)
List_Pointer(Problem_S.Group,
PostSubOperation_P->PostQuantitySupport[Order]))->InitialList ;
}
else {
PQ_P = NCPQ_P ; Support_L = NULL ;
}
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);
Group_P = (PostSubOperation_P->Case.OnRegion.RegionIndex < 0)? NULL :
(struct Group *)
List_Pointer(Problem_S.Group,
PostSubOperation_P->Case.OnRegion.RegionIndex);
Region_L = Group_P? Group_P->InitialList : NULL ;
Group_FunctionType = Group_P? Group_P->FunctionType : REGION;
if (!Support_L &&
List_Nbr(NCPQ_P->PostQuantityTerm) &&
(
((struct PostQuantityTerm *)List_Pointer(NCPQ_P->PostQuantityTerm, 0))
->Type == LOCALQUANTITY &&
((struct PostQuantityTerm *)List_Pointer(NCPQ_P->PostQuantityTerm, 0))
->EvaluationType == LOCAL)
) {
if (Group_FunctionType == REGION)
Msg(GERROR, "Print OnRegion not valid for PostProcessing Quantity '%s'",
NCPQ_P->Name);
else
Type_Evaluation = LOCAL;
}
else
Type_Evaluation = GLOBAL;
if (Region_L) {
if (Group_P->FunctionType == REGION) {
List_Sort(Region_L, fcmp_int) ;
Nbr_Region = List_Nbr(Region_L) ;
if (PostSubOperation_P->Format != FORMAT_SPACE_TABLE) {
fprintf(PostStream, "# %s on", PQ_P->Name) ;
for(i = 0 ; i < Nbr_Region ; i++) {
List_Read(Region_L, i, &Num_Region) ;
fprintf(PostStream, " %d", Num_Region) ;
}
fprintf(PostStream, "\n") ;
}
}
else if (Group_P->FunctionType == NODESOF) {
if (!Group_P->ExtendedList) Generate_ExtendedGroup(Group_P) ;
Region_L = Group_P->ExtendedList ; /* Attention: new Region_L */
Nbr_Region = List_Nbr(Region_L) ;
if (PostSubOperation_P->Comma) /* Provisoire */
Flag_Summation = 1;
}
else {
Msg(GERROR, "Function type (%d) not allowed for PrintOnRegion",
Group_P->FunctionType) ;
}
}
else
Nbr_Region = 1 ;
if (Type_Evaluation == LOCAL)
Init_SearchGrid(&Current.GeoData->Grid) ;
for (iTime = 0 ; iTime < NbrTimeStep ; iTime++) {
Pos_InitAllSolutions(PostSubOperation_P->TimeStep_L, iTime) ;
if (Flag_Summation) {
Cal_ZeroValue(&ValueSummed) ;
}
for(i = 0 ; i < Nbr_Region ; i++) {
if (Region_L)
List_Read(Region_L, i, &Num_Region) ;
else
Num_Region = NO_REGION ;
Current.SubRegion = Num_Region ; /* Region being a GlobalQuantity Entity no */
Current.NumEntity = Num_Region ; /* for OnRegion NodesOf */
Element.GeoElement = NULL ;
Element.Num = NO_ELEMENT ;
Element.Type = -1 ;
Current.Region = Element.Region = Num_Region ;
Current.x = Current.y = Current.z = 0. ;
if (Type_Evaluation == GLOBAL) {
Cal_PostQuantity(PQ_P, DefineQuantity_P0, QuantityStorage_P0,
Support_L, &Element, 0., 0., 0., &Value) ;
}
else {
if (Group_FunctionType == NODESOF)
Geo_GetNodesCoordinates(1, &Num_Region,
&Current.x, &Current.y, &Current.z) ;
InWhichElement(Current.GeoData->Grid, NULL, &Element,
PostSubOperation_P->Dimension,
Current.x, Current.y, Current.z, &u, &v, &w) ;
Cal_PostQuantity(PQ_P, DefineQuantity_P0, QuantityStorage_P0,
Support_L, &Element, u, v, w, &Value) ;
}
if (PostSubOperation_P->StoreInRegister >= 0)
Cal_StoreInRegister(&Value, PostSubOperation_P->StoreInRegister) ;
Format_PostValue(PostSubOperation_P->Format, PostSubOperation_P->Comma,
Group_FunctionType,
Current.Time, i, Current.NumEntity, Nbr_Region,
Current.NbrHar, PostSubOperation_P->HarmonicToTime,
PostSubOperation_P->NoNewLine,
&Value) ;
if (Flag_Summation) {
ValueSummed.Type = Value.Type ;
Cal_AddValue(&ValueSummed, &Value, &ValueSummed);
}
}
if (Flag_Summation && PostSubOperation_P->Format == FORMAT_REGION_TABLE) {
fprintf(PostStream, "#Sum: ") ;
Print_Value(&ValueSummed);
fprintf(PostStream, "\n") ;
}
}
Format_PostFooter(PostSubOperation_P, 0);
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* P o s _ P r i n t W i t h A r g u m e n t */
/* ------------------------------------------------------------------------ */
void Pos_PrintWithArgument(struct PostQuantity *NCPQ_P,
struct PostQuantity *CPQ_P,
int Order,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct PostSubOperation *PostSubOperation_P) {
struct Element Element ;
struct Value Value ;
struct Expression * Expression_P ;
List_T *Region_L ;
int i, N, Num_Region ;
double X[2], S, x ;
GetDP_Begin("Pos_PrintWithArgument");
if(CPQ_P)
Msg(GERROR, "Cumulative PostProcessing Quantity in PrintWithArgument not done") ;
X[0] = PostSubOperation_P->Case.WithArgument.x[0] ;
X[1] = PostSubOperation_P->Case.WithArgument.x[1] ;
N = PostSubOperation_P->Case.WithArgument.n ;
Expression_P = (struct Expression *)
List_Pointer(Problem_S.Expression,
PostSubOperation_P->Case.WithArgument.ArgumentIndex) ;
Region_L = ((struct Group *)
List_Pointer(Problem_S.Group,
PostSubOperation_P->Case.WithArgument.RegionIndex))
->InitialList ;
if (List_Nbr(Region_L))
List_Read(Region_L, 0, &Num_Region) ;
else
Num_Region = NO_REGION ;
for (i = 0 ; i <= N ; i++) {
S = (double)i / (double)(N ? N : 1) ;
x = X[0] + (X[1] - X[0]) * S ;
Expression_P->Case.Constant = x ;
Element.GeoElement = NULL ;
Element.Num = NO_ELEMENT ;
Element.Type = -1 ;
Current.Region = Element.Region = Num_Region ;
Current.x = Current.y = Current.z = 0. ;
Cal_PostQuantity(NCPQ_P, DefineQuantity_P0, QuantityStorage_P0,
NULL, &Element, 0., 0., 0., &Value) ;
Format_PostValue(PostSubOperation_P->Format, PostSubOperation_P->Comma,
REGION,
x, 0, 0, 1,
Current.NbrHar, PostSubOperation_P->HarmonicToTime,
PostSubOperation_P->NoNewLine,
&Value) ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* P o s _ P r i n t E x p r e s s i o n */
/* ------------------------------------------------------------------------ */
void Pos_PrintExpression(struct PostSubOperation *PostSubOperation_P){
int NbrTimeStep, iTime;
struct Value Value;
char *str = PostSubOperation_P->Case.Expression.String;
char *str2 = PostSubOperation_P->Case.Expression.String2;
int expr = PostSubOperation_P->Case.Expression.ExpressionIndex;
GetDP_Begin("Pos_PrintExpression");
NbrTimeStep = Pos_InitTimeSteps(PostSubOperation_P);
for(iTime = 0; iTime < NbrTimeStep; iTime++){
Pos_InitAllSolutions(PostSubOperation_P->TimeStep_L, iTime) ;
if(expr >= 0){
Get_ValueOfExpressionByIndex(expr, NULL, 0., 0., 0., &Value) ;
if(str) fprintf(PostStream, "%s", str);
Print_Value(&Value);
fprintf(PostStream, "\n") ;
}
else if(str2){
if(str) fprintf(PostStream, "%s", str);
fprintf(PostStream, "%s\n", str2);
}
else if(str){
fprintf(PostStream, "%s\n", str);
}
}
GetDP_End;
}
/* ------------------------------------------------------------------------ */
/* P o s _ P r i n t G r o u p */
/* ------------------------------------------------------------------------ */
void Pos_PrintGroup(struct PostSubOperation *PostSubOperation_P) {
struct Group *Group_P;
struct Geo_Element *GeoElement;
struct PostElement *SL;
List_T *Region_L;
int i, NbrGeo, iGeo, *NumNodes;
double x [NBR_MAX_NODES_IN_ELEMENT] ;
double y [NBR_MAX_NODES_IN_ELEMENT] ;
double z [NBR_MAX_NODES_IN_ELEMENT] ;
GetDP_Begin("Pos_PrintGroup");
NbrGeo = Geo_GetNbrGeoElements() ;
Format_PostHeader(PostSubOperation_P->Format,
PostSubOperation_P->Iso, 1,
PostSubOperation_P->HarmonicToTime,
PostSubOperation_P->CombinationType, 0,
NULL, NULL);
Region_L = ((struct Group *)
List_Pointer(Problem_S.Group,
PostSubOperation_P->Case.Group.GroupIndex))->InitialList ;
Group_P = (struct Group *)
List_Pointer(Problem_S.Group,
PostSubOperation_P->Case.Group.ExtendedGroupIndex);
SL = Create_PostElement(0, LINE, 2, 1) ;
if(!Group_P->ExtendedList) Generate_ExtendedGroup(Group_P) ;
for(iGeo = 0 ; iGeo < NbrGeo ; iGeo++) {
if(InteractiveInterrupt) break;
Progress(iGeo, NbrGeo, "Compute: ") ;
GeoElement = Geo_GetGeoElement(iGeo) ;
if(List_Search(Region_L, &GeoElement->Region, fcmp_int)){
Geo_GetNodesCoordinates
(GeoElement->NbrNodes, GeoElement->NumNodes, x, y, z) ;
switch (Group_P->FunctionType) {
case EDGESOFTREEIN :
if(!GeoElement->NbrEdges) Geo_CreateEdgesOfElement(GeoElement) ;
for(i=0 ; i<GeoElement->NbrEdges ; i++){
if(List_Search(Group_P->ExtendedList, &GeoElement->NumEdges[i], fcmp_absint)){
NumNodes = Geo_GetNodesOfEdgeInElement(GeoElement, i) ;
SL->Index = iGeo;
SL->x[0] = x[abs(NumNodes[0])-1]; SL->x[1] = x[abs(NumNodes[1])-1];
SL->y[0] = y[abs(NumNodes[0])-1]; SL->y[1] = y[abs(NumNodes[1])-1];
SL->z[0] = z[abs(NumNodes[0])-1]; SL->z[1] = z[abs(NumNodes[1])-1];
SL->Value[0].Type = SL->Value[1].Type = SCALAR ;
SL->Value[0].Val[0] = SL->Value[1].Val[0] = GeoElement->NumEdges[i];
Format_PostElement(PostSubOperation_P->Format, PostSubOperation_P->Iso, 0,
0, 0, 1, 1, 1,
NULL, SL,
PostSubOperation_P->ChangeOfCoordinates,
PostSubOperation_P->ChangeOfValues);
}
}
break ;
default :
Msg(GERROR, "Print function not implemented for this kind of Group");
break ;
}
}
}
Destroy_PostElement(SL) ;
Format_PostFooter(PostSubOperation_P, 0);
GetDP_End ;
}
syntax highlighted by Code2HTML, v. 0.9.1