#define RCSID "$Id: Pos_Quantity.c,v 1.20 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 "Treatment_Formulation.h"
#include "Pos_Formulation.h"
#include "Pos_Quantity.h"
#include "Get_DofOfElement.h"
#include "GeoData.h"
#include "Cal_Quantity.h"
#include "Get_Geometry.h"
#include "Data_Passive.h"
#include "Data_DefineE.h"
#include "CurrentData.h"
#include "Tools.h"
/* ------------------------------------------------------------------------ */
/* C a l _ P o s t Q u a n t i t y */
/* ------------------------------------------------------------------------ */
void Cal_PostQuantity(struct PostQuantity *PostQuantity_P,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
List_T *Support_L,
struct Element *Element,
double u, double v, double w,
struct Value *Value) {
struct PostQuantityTerm PostQuantityTerm ;
List_T *InRegion_L ;
int i_PQT, Type_Quantity ;
GetDP_Begin("Cal_PostQuantity");
/* mettre tout a zero: on ne connait pas a priori le type de retour */
/* (default type and value returned if Type_Quantity == -1) */
Cal_ZeroValue(Value);
Value->Type = SCALAR;
/* Loop on PostQuantity Terms */
/* ... with sum of results if common supports (In ...) */
for (i_PQT = 0 ; i_PQT < List_Nbr(PostQuantity_P->PostQuantityTerm) ; i_PQT++) {
List_Read(PostQuantity_P->PostQuantityTerm, i_PQT, &PostQuantityTerm) ;
InRegion_L = (PostQuantityTerm.InIndex < 0)? NULL :
((struct Group *)List_Pointer(Problem_S.Group,
PostQuantityTerm.InIndex))->InitialList ;
if (!Support_L) Type_Quantity = PostQuantityTerm.Type ;
else Type_Quantity = GLOBALQUANTITY ; /* Always if Support */
if (InRegion_L) {
if (Element->Num != NO_ELEMENT) {
if (!List_Search(InRegion_L, &Element->Region, fcmp_int)) {
Type_Quantity = -1 ;
}
}
else {
if (Type_Quantity == GLOBALQUANTITY) {
/* Plus de test ici... vu que le OnRegion de la PostOperation n'a rien
a voir avec le support d'une integration ...
if (!List_Search(InRegion_L, &Current.Region, fcmp_int)) {
Type_Quantity = -1 ;
}
*/
/* Il faut plutot voir si il existe au moins une region de InRegion_L
qui soit dans Support_L ... cela est fait apres, pour chaque element */
}
else if (Type_Quantity != INTEGRALQUANTITY) {
Type_Quantity = -1 ;
}
}
}
/* else if !InRegion_L -> No filter, i.e. globally defined quantity */
/* ---------------------------- */
/* Local or Integral quantities */
/* ---------------------------- */
if (Type_Quantity == LOCALQUANTITY || Type_Quantity == INTEGRALQUANTITY) {
Pos_LocalOrIntegralQuantity(PostQuantity_P,
DefineQuantity_P0, QuantityStorage_P0,
&PostQuantityTerm, Element, Type_Quantity,
u, v, w, Value) ;
}
/* ----------------- */
/* Global quantities */
/* ----------------- */
else if (Type_Quantity == GLOBALQUANTITY) {
Pos_GlobalQuantity(PostQuantity_P,
DefineQuantity_P0, QuantityStorage_P0,
&PostQuantityTerm, Element, InRegion_L, Support_L, Value) ;
}
} /* for i_PQT ... */
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* P o s _ G l o b a l Q u a n t i t y */
/* ------------------------------------------------------------------------ */
void Pos_GlobalQuantity(struct PostQuantity *PostQuantity_P,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct PostQuantityTerm *PostQuantityTerm_P,
struct Element *ElementEmpty,
List_T *InRegion_L, List_T *Support_L,
struct Value *Value) {
struct DefineQuantity *DefineQuantity_P ;
struct QuantityStorage *QuantityStorage_P ;
struct FunctionSpace *FunctionSpace_P ;
struct GlobalQuantity *GlobalQuantity_P ;
struct Value TermValue ;
int k, Index_DefineQuantity ;
int Nbr_Element, i_Element ;
struct Element Element ;
int Type_Quantity ;
GetDP_Begin("Pos_GlobalQuantity");
if (PostQuantityTerm_P->EvaluationType == LOCAL &&
List_Search(InRegion_L, &Current.Region, fcmp_int)) {
for (k = 0 ; k < PostQuantityTerm_P->NbrQuantityIndex ; k++) {
Index_DefineQuantity = PostQuantityTerm_P->QuantityIndexTable[k] ;
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
if (QuantityStorage_P->NumLastElementForFunctionSpace != Current.Region) {
QuantityStorage_P->NumLastElementForFunctionSpace = Current.Region ;
QuantityStorage_P->FunctionSpace = FunctionSpace_P =
(struct FunctionSpace*)
List_Pointer(Problem_S.FunctionSpace,
DefineQuantity_P->FunctionSpaceIndex) ;
GlobalQuantity_P = (struct GlobalQuantity*)
List_Pointer
(QuantityStorage_P->FunctionSpace->GlobalQuantity,
*(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
if (DefineQuantity_P->Type == GLOBALQUANTITY)
Get_DofOfRegion(Current.Region, GlobalQuantity_P,
FunctionSpace_P, QuantityStorage_P) ;
}
}
Cal_WholeQuantity
(Current.Element = ElementEmpty,
QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
Current.u = 0., Current.v = 0., Current.w = 0., -1, -1, &TermValue) ;
Value->Type = TermValue.Type;
Cal_AddValue(Value,&TermValue,Value);
} /* if LOCAL && ... */
else if (PostQuantityTerm_P->EvaluationType == INTEGRAL) {
Nbr_Element = Geo_GetNbrGeoElements() ;
Get_InitDofOfElement(&Element) ;
Type_Quantity = LOCALQUANTITY ; /* Attention... il faut se comprendre: */
/* il s'agit de grandeurs locales qui seront integrees */
for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
Progress(i_Element, Nbr_Element, "Accumulate: ") ;
Element.GeoElement = Geo_GetGeoElement(i_Element) ;
Element.Num = Element.GeoElement->Num ;
Element.Type = Element.GeoElement->Type ;
Current.Region = Element.Region = Element.GeoElement->Region ;
/* Filter: only elements in both InRegion_L and Support_L are considered */
if ((!InRegion_L || List_Search(InRegion_L, &Element.Region, fcmp_int)) &&
(!Support_L || List_Search(Support_L, &Element.Region, fcmp_int))) {
Get_NodesCoordinatesOfElement(&Element) ;
Current.x = Element.x[0];
Current.y = Element.y[0];
Current.z = Element.z[0];
Pos_LocalOrIntegralQuantity(PostQuantity_P,
DefineQuantity_P0, QuantityStorage_P0,
PostQuantityTerm_P, &Element, Type_Quantity,
0., 0., 0., Value) ;
}
} /* for i_Element ... */
} /* if INTEGRAL ... */
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* C a l _ P o s t C u m u l a t i v e Q u a n t i t y */
/* ------------------------------------------------------------------------ */
void Cal_PostCumulativeQuantity(List_T *Region_L,
int SupportIndex,
List_T *TimeStep_L,
struct PostQuantity *PostQuantity_P,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct Value **Values) {
struct Element Element ;
List_T *Support_L ;
int i, NbrTimeStep ;
GetDP_Begin("Cal_PostCumulativeQuantity");
Support_L = ((struct Group *)
List_Pointer(Problem_S.Group, SupportIndex))->InitialList ;
NbrTimeStep = List_Nbr(TimeStep_L) ;
*Values = (struct Value *)Malloc(NbrTimeStep*sizeof(struct Value)) ;
Element.Num = NO_ELEMENT ;
for(i = 0 ; i < NbrTimeStep ; i++) {
Pos_InitAllSolutions(TimeStep_L, i) ;
Cal_PostQuantity(PostQuantity_P, DefineQuantity_P0, QuantityStorage_P0,
Support_L, &Element, 0, 0, 0, &(*Values)[i]) ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* C o m b i n e _ P o s t Q u a n t i t y */
/* ------------------------------------------------------------------------ */
void Combine_PostQuantity(int Type, int Order,
struct Value *V1, struct Value *V2){
GetDP_Begin("Combine_PostQuantity");
switch(Type){
case MULTIPLICATION : Cal_ProductValue(V1, V2, V1) ; break ;
case ADDITION : Cal_AddValue(V1, V2, V1) ; break ;
case DIVISION : Cal_DivideValue(Order?V1:V2, Order?V2:V1, V1) ; break;
case SOUSTRACTION : Cal_SubstractValue(Order?V1:V2, Order?V2:V1, V1) ; break;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* P o s _ L o c a l O r I n t e g r a l Q u a n t i t y */
/* ------------------------------------------------------------------------ */
static int Warning_NoJacobian = 0 ;
void Pos_LocalOrIntegralQuantity(struct PostQuantity *PostQuantity_P,
struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage *QuantityStorage_P0,
struct PostQuantityTerm *PostQuantityTerm_P,
struct Element *Element,
int Type_Quantity,
double u, double v, double w,
struct Value *Value) {
struct FunctionSpace *FunctionSpace_P ;
struct DefineQuantity *DefineQuantity_P ;
struct QuantityStorage *QuantityStorage_P ;
struct Value TermValue, tmpValue;
struct JacobianCase *JacobianCase_P0 ;
struct IntegrationCase *IntegrationCase_P ;
struct Quadrature *Quadrature_P ;
List_T *IntegrationCase_L, *JacobianCase_L ;
double ui, vi, wi, weight, Factor ;
int Index_DefineQuantity ;
int i, j, Type_Dimension ;
int CriterionIndex, Nbr_IntPoints, i_IntPoint ;
double (*Get_Jacobian) (struct Element * Element, MATRIX3x3 * Jac) = 0;
void (*Get_IntPoint) (int Nbr_Points, int Num,
double * u, double * v, double * w, double * wght) ;
GetDP_Begin("Pos_LocalOrIntegralQuantity");
/* Get the functionspaces
Get the DoF for local quantities */
for (i = 0 ; i < PostQuantityTerm_P->NbrQuantityIndex ; i++) {
Index_DefineQuantity = PostQuantityTerm_P->QuantityIndexTable[i] ;
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
if (QuantityStorage_P->NumLastElementForFunctionSpace != Element->Num) {
QuantityStorage_P->NumLastElementForFunctionSpace = Element->Num ;
if (Type_Quantity != INTEGRALQUANTITY){
QuantityStorage_P->FunctionSpace = FunctionSpace_P =
(struct FunctionSpace*)
List_Pointer(Problem_S.FunctionSpace,
DefineQuantity_P->FunctionSpaceIndex) ;
if (DefineQuantity_P->Type == LOCALQUANTITY)
Get_DofOfElement(Element, FunctionSpace_P, QuantityStorage_P,
DefineQuantity_P->IndexInFunctionSpace) ;
}
else{ /* INTEGRALQUANTITY */
if(DefineQuantity_P->IntegralQuantity.DefineQuantityIndexDof >= 0)
QuantityStorage_P->FunctionSpace = (struct FunctionSpace*)
List_Pointer(Problem_S.FunctionSpace,
DefineQuantity_P->FunctionSpaceIndex) ;
/* Get the function space for the associated local quantities */
for (j = 0 ; j < DefineQuantity_P->IntegralQuantity.NbrQuantityIndex ; j++) {
Index_DefineQuantity = DefineQuantity_P->IntegralQuantity.QuantityIndexTable[j];
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
QuantityStorage_P->FunctionSpace = (struct FunctionSpace*)
List_Pointer(Problem_S.FunctionSpace,
DefineQuantity_P->FunctionSpaceIndex) ;
}
}
}
}
/* get the jacobian */
if (Element->Num != NO_ELEMENT) {
if (PostQuantityTerm_P->JacobianMethodIndex >= 0) {
JacobianCase_L =
((struct JacobianMethod *)
List_Pointer(Problem_S.JacobianMethod,
PostQuantityTerm_P->JacobianMethodIndex))
->JacobianCase ;
JacobianCase_P0 = (struct JacobianCase*)List_Pointer(JacobianCase_L, 0);
i = 0 ;
while ((i < List_Nbr(JacobianCase_L)) &&
((j = (JacobianCase_P0 + i)->RegionIndex) >= 0) &&
!List_Search
(((struct Group *)List_Pointer(Problem_S.Group, j))
->InitialList, &Element->Region, fcmp_int) ) i++ ;
if (i == List_Nbr(JacobianCase_L))
Msg(GERROR, "Undefined Jacobian in Region %d", Element->Region) ;
Element->JacobianCase = JacobianCase_P0 + i ;
Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
Element->Type, &Type_Dimension) ;
}
else {
if(!Warning_NoJacobian){
Msg(WARNING, "No Jacobian method specification in PostProcessing quantity");
Msg(WARNING, "Using default Jacobian (Vol)");
Warning_NoJacobian = 1 ;
}
Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
Get_JacobianFunction(JACOBIAN_VOL, Element->Type, &Type_Dimension) ;
}
Get_NodesCoordinatesOfElement(Element) ;
}
/* local evaluation at one point */
if(PostQuantityTerm_P->EvaluationType == LOCAL){
if (Element->Num != NO_ELEMENT) {
Get_BFGeoElement(Element, u, v, w) ;
Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
if (Element->DetJac != 0.)
Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
&Element->Jac, &Element->InvJac) ;
}
Cal_WholeQuantity
(Current.Element = Element,
QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
Current.u = u, Current.v = v, Current.w = w, -1, -1, &TermValue) ;
}
/* integral evaluation over the element */
else if(PostQuantityTerm_P->EvaluationType == INTEGRAL){
if(Element->Num == NO_ELEMENT)
Msg(GERROR, "No element in which to integrate");
if(PostQuantityTerm_P->IntegrationMethodIndex < 0)
Msg(GERROR, "Missing Integration method in PostProcesssing Quantity '%s'",
PostQuantity_P->Name);
IntegrationCase_L =
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
PostQuantityTerm_P->IntegrationMethodIndex))->IntegrationCase ;
CriterionIndex =
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
PostQuantityTerm_P->IntegrationMethodIndex))
->CriterionIndex ;
IntegrationCase_P = Get_IntegrationCase(Element,
IntegrationCase_L,
CriterionIndex) ;
if(IntegrationCase_P->Type != GAUSS)
Msg(GERROR, "Only numerical integration is available "
"in Integral PostQuantities");
Quadrature_P = (struct Quadrature*)
List_PQuery(IntegrationCase_P->Case, &Element->Type, fcmp_int);
if(!Quadrature_P)
Msg(GERROR, "Unknown type of Element (%s) for Integration method (%s) "
" in PostProcessing Quantity (%s)",
Get_StringForDefine(Element_Type, Element->Type),
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
PostQuantityTerm_P->IntegrationMethodIndex))->Name,
PostQuantity_P->Name);
Cal_ZeroValue(&TermValue);
Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
Get_IntPoint = (void (*) (int,int,double*,double*,double*,double*))
Quadrature_P->Function ;
for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
Get_IntPoint(Nbr_IntPoints, i_IntPoint, &ui, &vi, &wi, &weight) ;
Get_BFGeoElement (Element, ui, vi, wi) ;
Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
if (Element->DetJac != 0.){
Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
&Element->Jac, &Element->InvJac) ;
}
else{
Msg(WARNING, "Zero determinant in 'Cal_PostQuantity'");
}
Current.x = Current.y = Current.z = 0. ;
if (Type_Quantity == INTEGRALQUANTITY){
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
Current.x += Element->x[i] * Element->n[i] ;
Current.y += Element->y[i] * Element->n[i] ;
Current.z += Element->z[i] * Element->n[i] ;
}
}
Cal_WholeQuantity
(Current.Element = Element,
QuantityStorage_P0, PostQuantityTerm_P->WholeQuantity,
Current.u = ui, Current.v = vi, Current.w = wi, -1, -1, &tmpValue) ;
Factor = weight * fabs(Element->DetJac) ;
TermValue.Type = tmpValue.Type ;
Cal_AddMultValue(&TermValue,&tmpValue,Factor,&TermValue);
}
}
Value->Type = TermValue.Type;
Cal_AddValue(Value,&TermValue,Value);
GetDP_End ;
}
syntax highlighted by Code2HTML, v. 0.9.1