#define RCSID "$Id: Pos_FemInterpolation.c,v 1.21 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 . */ #include "GetDP.h" #include "Treatment_Formulation.h" #include "Get_DofOfElement.h" #include "Get_Geometry.h" #include "GeoData.h" #include "DofData.h" #include "Pos_Search.h" #include "CurrentData.h" #include "Tools.h" /* ------------------------------------------------------------------------ */ /* P o s _ F e m I n t e r p o l a t i o n */ /* ------------------------------------------------------------------------ */ int SearchType(int type){ GetDP_Begin("SearchType"); switch(type){ case POINT : GetDP_Return(_ALL) ; case LINE : GetDP_Return(_1D) ; case TRIANGLE : case QUADRANGLE : GetDP_Return(_2D) ; case TETRAHEDRON : case HEXAHEDRON : case PRISM : case PYRAMID : GetDP_Return(_3D) ; default : Msg(GERROR, "Unknown Type in 'SearchType'"); GetDP_Return(-1) ; } } void Pos_FemInterpolation(struct Element * Element, struct QuantityStorage * QuantityStorage_P0, struct QuantityStorage * QuantityStorage_P, int Type_Quantity, int Type_Operator, int Type_Dimension, int UseXYZ, double u, double v, double w, double x, double y, double z, double Val[], int * Type_Value, int Flag_ChangeOfCoordinates) { void (*xFunctionBF[NBR_MAX_BASISFUNCTIONS]) (struct Element *, int, double, double, double, double []) ; void (*xChangeOfCoordinates) () = 0; struct IntegralQuantityActive IQA ; struct Value vBFxDof[NBR_MAX_BASISFUNCTIONS] ; struct GeoData * GeoData_P ; struct Element TheElement, * TheElement_P ; struct QuantityStorage * QS_P ; double vBFu[NBR_MAX_BASISFUNCTIONS][MAX_DIM] ; double Val_Dof, Val_Dof_r, Val_Dof_i ; int Type_DefineQuantity, SubType_DefineQuantity, Type_Form ; int i, j, k, Nbr_Dof = 0 ; int GeoDataNum = 0, UseNewGeo = 0 ; GetDP_Begin("Pos_FemInterpolation"); /* ------------- Quantity Type ------------- */ Type_DefineQuantity = QuantityStorage_P->DefineQuantity->Type ; if(Type_DefineQuantity == INTEGRALQUANTITY){ if(QuantityStorage_P->DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < 0){ SubType_DefineQuantity = NODOF ; } else{ SubType_DefineQuantity = INTEGRALQUANTITY ; } } else{ SubType_DefineQuantity = Type_DefineQuantity ; } /* --------------- Get The Element --------------- */ if(SubType_DefineQuantity != NODOF) { if(!QuantityStorage_P->FunctionSpace) Msg(GERROR, "No available function space for quantity"); if(!QuantityStorage_P->FunctionSpace->DofData) Msg(GERROR, "No available data to interpolate quantity"); GeoDataNum = QuantityStorage_P->FunctionSpace->DofData->GeoDataIndex; UseNewGeo = (GeoDataNum != Current.GeoData->Num) ; if(UseXYZ || UseNewGeo){ if(UseNewGeo){ GeoData_P = (struct GeoData *)List_Pointer(GeoData_L, GeoDataNum); GeoDataNum = Current.GeoData->Num ; Geo_SetCurrentGeoData(Current.GeoData = GeoData_P) ; } if(!UseXYZ){ x = y = z = 0. ; for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) { x += Element->x[i] * Element->n[i] ; y += Element->y[i] * Element->n[i] ; z += Element->z[i] * Element->n[i] ; } } Init_SearchGrid(&Current.GeoData->Grid); InWhichElement(Current.GeoData->Grid, NULL, &TheElement, (Type_Dimension >= 0) ? Type_Dimension : SearchType(Element->Type), x, y, z, &u, &v, &w) ; TheElement_P = &TheElement ; Get_InitDofOfElement(&TheElement) ; Get_DofOfElement (&TheElement, QuantityStorage_P->FunctionSpace, QuantityStorage_P, QuantityStorage_P->DefineQuantity->IndexInFunctionSpace) ; } else{ TheElement_P = Element; } } else{ TheElement_P = Element ; } /* ------------------ Init LocalQuantity ------------------ */ if (Type_DefineQuantity == LOCALQUANTITY) { if (TheElement_P->Num != NO_ELEMENT) { Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction ; Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBF, Type_Operator, QuantityStorage_P, &Type_Form) ; xChangeOfCoordinates = (void (*)())Get_ChangeOfCoordinates ((Flag_ChangeOfCoordinates && TheElement_P->Num != NO_ELEMENT), Type_Form) ; } else { Msg(WARNING, "No element found in mesh for LocalQuantity interpolation"); Nbr_Dof = 0 ; Type_Form = VECTOR ; } } /* --------------------- Init IntegralQuantity --------------------- */ else if (Type_DefineQuantity == INTEGRALQUANTITY) { if(Type_Operator != NOOP){ Msg(GERROR, "Operator acting on Integral Quantity"); } Type_Form = VECTOR ; Get_InitElementSource(TheElement_P, QuantityStorage_P->DefineQuantity->IntegralQuantity.InIndex) ; IQA.IntegrationCase_L = ((struct IntegrationMethod *) List_Pointer(Problem_S.IntegrationMethod, QuantityStorage_P->DefineQuantity-> IntegralQuantity.IntegrationMethodIndex)) ->IntegrationCase ; IQA.CriterionIndex = ((struct IntegrationMethod *) List_Pointer(Problem_S.IntegrationMethod, QuantityStorage_P->DefineQuantity-> IntegralQuantity.IntegrationMethodIndex)) ->CriterionIndex ; IQA.JacobianCase_L = ((struct JacobianMethod *) List_Pointer(Problem_S.JacobianMethod, QuantityStorage_P->DefineQuantity-> IntegralQuantity.JacobianMethodIndex)) ->JacobianCase ; xChangeOfCoordinates = (void (*)())Get_ChangeOfCoordinates(0, Type_Form) ; } /* ---------------------- Compute GlobalQuantity ---------------------- */ if (Type_DefineQuantity == GLOBALQUANTITY) { if(Current.NbrHar==1){ if (Type_Quantity == QUANTITY_BF) Val[0] = (QuantityStorage_P->BasisFunction[0].Dof->Entity == Current.SubRegion)? 1. : 0. ; else Dof_GetRealDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[0].Dof, &Val[0]) ; } else{ for (k = 0 ; k < Current.NbrHar ; k+=2) { if (Type_Quantity == QUANTITY_BF) { Val[MAX_DIM*k] = (QuantityStorage_P->BasisFunction[0].Dof->Entity == Current.SubRegion)? 1. : 0. ; Val[MAX_DIM*(k+1)] = 0. ; } else { Dof_GetComplexDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[0].Dof + k/2*gCOMPLEX_INCREMENT, &Val[MAX_DIM*k], &Val[MAX_DIM*(k+1)]) ; } } } *Type_Value = SCALAR ; GetDP_End ; } /* ----------------------------------- Compute Local / Integral Quantities ----------------------------------- */ i = Current.NbrHar * MAX_DIM ; for (k = 0 ; k < i ; k++) Val[k] = 0. ; while (1) { if (Type_DefineQuantity == INTEGRALQUANTITY) { if (Get_NextElementSource(TheElement_P->ElementSource)) { Get_NodesCoordinatesOfElement(TheElement_P->ElementSource) ; if(SubType_DefineQuantity != NODOF){ Get_DofOfElement(TheElement_P->ElementSource, QuantityStorage_P->FunctionSpace, QuantityStorage_P, QuantityStorage_P->DefineQuantity->IndexInFunctionSpace) ; Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction ; Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBF, QuantityStorage_P->DefineQuantity->IntegralQuantity.TypeOperatorDof, QuantityStorage_P, &IQA.Type_FormDof) ; Type_Form = IQA.Type_FormDof ; /* good form */ } else{ Nbr_Dof = 1 ; xFunctionBF[0] = NULL ; /* for analytic integration tests */ Type_Form = IQA.Type_FormDof = VECTOR ; /* form type unknown */ for (j = 0 ; j < QuantityStorage_P->DefineQuantity->IntegralQuantity.NbrQuantityIndex ; j++) { QS_P = QuantityStorage_P0 + QuantityStorage_P->DefineQuantity->IntegralQuantity.QuantityIndexTable[j] ; Get_DofOfElement(TheElement_P->ElementSource, QS_P->FunctionSpace, QS_P, QS_P->DefineQuantity->IndexInFunctionSpace) ; } } Cal_InitIntegralQuantity (TheElement_P, &IQA, QuantityStorage_P); } else break ; } /* ----- Local ----- */ if (Type_DefineQuantity == LOCALQUANTITY) { if (TheElement_P->Num != NO_ELEMENT) { for (j = 0 ; j < Nbr_Dof ; j++) { xFunctionBF[j] (TheElement_P, QuantityStorage_P->BasisFunction[j].NumEntityInElement+1, u, v, w, vBFu[j]) ; ((void (*)(struct Element*, double*, double*)) xChangeOfCoordinates) (TheElement_P, vBFu[j], vBFxDof[j].Val) ; /* printf("j %d , Num %d, Type_Form %d change %d " "vBFu[j] %e %e %e vBFx[j] %e %e %e\n", j, QuantityStorage_P->BasisFunction[j].NumEntityInElement+1, Type_Form,(Flag_ChangeOfCoordinates && TheElement_P->Num != NO_ELEMENT), vBFu[j][0],vBFu[j][1],vBFu[j][2], vBFxDof[j].Val[0],vBFxDof[j].Val[1],vBFxDof[j].Val[2]); */ } } /* interpolate (vBFxDof is real-valued) */ switch (Type_Form) { case FORM0 : case FORM3 : case FORM3P : case SCALAR : if(Current.NbrHar==1){ for (j = 0 ; j < Nbr_Dof ; j++){ if (Type_Quantity == QUANTITY_BF) Val_Dof = (QuantityStorage_P->BasisFunction[j].Dof->Entity == Current.SubRegion)? 1. : 0. ; else Dof_GetRealDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof, &Val_Dof) ; Val[0] += vBFxDof[j].Val[0] * Val_Dof ; } } else{ for (j = 0 ; j < Nbr_Dof ; j++){ for (k = 0 ; k < Current.NbrHar ; k+=2) { if (Type_Quantity == QUANTITY_BF) { Val_Dof_r = (QuantityStorage_P->BasisFunction[j].Dof->Entity == Current.SubRegion)? 1. : 0. ; Val_Dof_i = 0. ; } else { Dof_GetComplexDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT, &Val_Dof_r, &Val_Dof_i) ; } Val[MAX_DIM*k] += vBFxDof[j].Val[0] * Val_Dof_r ; Val[MAX_DIM*(k+1)] += vBFxDof[j].Val[0] * Val_Dof_i ; } } } *Type_Value = SCALAR ; break ; case FORM1 : case FORM1P : case FORM2 : case FORM2P : case FORM1S : case FORM2S : case VECTOR : case VECTORP : if(Current.NbrHar==1){ for (j = 0 ; j < Nbr_Dof ; j++){ if (Type_Quantity == QUANTITY_BF) Val_Dof = (QuantityStorage_P->BasisFunction[j].Dof->Entity == Current.SubRegion)? 1. : 0. ; else Dof_GetRealDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof, &Val_Dof) ; Val[0] += vBFxDof[j].Val[0] * Val_Dof ; Val[1] += vBFxDof[j].Val[1] * Val_Dof ; Val[2] += vBFxDof[j].Val[2] * Val_Dof ; } } else{ for (j = 0 ; j < Nbr_Dof ; j++){ for (k = 0 ; k < Current.NbrHar ; k+=2) { if (Type_Quantity == QUANTITY_BF) { Val_Dof_r = (QuantityStorage_P->BasisFunction[j].Dof->Entity == Current.SubRegion)? 1. : 0. ; Val_Dof_i = 0. ; } else { Dof_GetComplexDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT, &Val_Dof_r, &Val_Dof_i) ; } Val[MAX_DIM*k ] += vBFxDof[j].Val[0] * Val_Dof_r ; Val[MAX_DIM*k+1] += vBFxDof[j].Val[1] * Val_Dof_r ; Val[MAX_DIM*k+2] += vBFxDof[j].Val[2] * Val_Dof_r ; Val[MAX_DIM*(k+1) ] += vBFxDof[j].Val[0] * Val_Dof_i ; Val[MAX_DIM*(k+1)+1] += vBFxDof[j].Val[1] * Val_Dof_i ; Val[MAX_DIM*(k+1)+2] += vBFxDof[j].Val[2] * Val_Dof_i ; } } } *Type_Value = VECTOR ; break ; default : Msg(GERROR, "Unknown Form type in 'Pos_FemInterpolation'"); break; } } /* -------- Integral -------- */ /* Ceci, c'est nul a chier. Ce qu'il faut faire, c'est ne pas reinterpoler ici, mais laisser au Cal_Quantity dans Cal_IntegralQuantity le soin de reinterpoler directment la quantity local intervenant ds la qte integrale s'il y a lieu ! Mais, comment faire avec l'integration analytique ? */ else { if (IQA.IntegrationCase_P->Type == ANALYTIC) Cal_AnalyticIntegralQuantity (Current.Element = TheElement_P, QuantityStorage_P, Nbr_Dof, (void (**)())xFunctionBF, vBFxDof) ; else Cal_NumericalIntegralQuantity (Current.Element = TheElement_P, &IQA, QuantityStorage_P0, QuantityStorage_P, SubType_DefineQuantity, Nbr_Dof, (void (**)())xFunctionBF, vBFxDof) ; Type_Form = vBFxDof[0].Type ; /* interpolate (vBFxDof can be complex-valued) */ if(SubType_DefineQuantity == NODOF){ switch (Type_Form) { case FORM0 : case FORM3 : case FORM3P : case SCALAR : for (k = 0 ; k < Current.NbrHar ; k++) Val[MAX_DIM*k] += vBFxDof[0].Val[MAX_DIM*k] ; *Type_Value = SCALAR ; break ; case FORM1 : case FORM1P : case FORM2 : case FORM2P : case FORM1S : case FORM2S : case VECTOR : case VECTORP : for (k = 0 ; k < Current.NbrHar ; k++) { Val[MAX_DIM*k] += vBFxDof[0].Val[MAX_DIM*k] ; Val[MAX_DIM*k+1] += vBFxDof[0].Val[MAX_DIM*k+1] ; Val[MAX_DIM*k+2] += vBFxDof[0].Val[MAX_DIM*k+2] ; } *Type_Value = VECTOR ; break ; default : Msg(GERROR, "Unknown Form type in 'Pos_FemInterpolation'"); break; } } else{ switch (Type_Form) { case FORM0 : case FORM3 : case FORM3P : case SCALAR : if(Current.NbrHar==1){ for (j = 0 ; j < Nbr_Dof ; j++){ Dof_GetRealDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof, &Val_Dof) ; Val[0] += vBFxDof[j].Val[0] * Val_Dof ; } } else{ for (j = 0 ; j < Nbr_Dof ; j++){ for (k = 0 ; k < Current.NbrHar ; k+=2) { Dof_GetComplexDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT, &Val_Dof_r, &Val_Dof_i) ; Val[MAX_DIM*k] += vBFxDof[j].Val[MAX_DIM*k] * Val_Dof_r - vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_i ; Val[MAX_DIM*(k+1)] += vBFxDof[j].Val[MAX_DIM*k] * Val_Dof_i + vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_r ; } } } *Type_Value = SCALAR ; break ; case FORM1 : case FORM1P : case FORM2 : case FORM2P : case FORM1S : case FORM2S : case VECTOR : case VECTORP : if(Current.NbrHar==1){ for (j = 0 ; j < Nbr_Dof ; j++){ Dof_GetRealDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof, &Val_Dof) ; Val[0] += vBFxDof[j].Val[0] * Val_Dof ; Val[1] += vBFxDof[j].Val[1] * Val_Dof ; Val[2] += vBFxDof[j].Val[2] * Val_Dof ; } } else{ for (j = 0 ; j < Nbr_Dof ; j++){ for (k = 0 ; k < Current.NbrHar ; k+=2) { Dof_GetComplexDofValue (QuantityStorage_P->FunctionSpace->DofData, QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT, &Val_Dof_r, &Val_Dof_i) ; Val[MAX_DIM*k] += vBFxDof[j].Val[MAX_DIM*k] * Val_Dof_r - vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_i ; Val[MAX_DIM*(k+1)] += vBFxDof[j].Val[MAX_DIM*k] * Val_Dof_i + vBFxDof[j].Val[MAX_DIM*(k+1)] * Val_Dof_r ; Val[MAX_DIM*k+1] += vBFxDof[j].Val[MAX_DIM*k+1] * Val_Dof_r - vBFxDof[j].Val[MAX_DIM*(k+1)+1] * Val_Dof_i ; Val[MAX_DIM*(k+1)+1] += vBFxDof[j].Val[MAX_DIM*k+1] * Val_Dof_i + vBFxDof[j].Val[MAX_DIM*(k+1)+1] * Val_Dof_r ; Val[MAX_DIM*k+2] += vBFxDof[j].Val[MAX_DIM*k+2] * Val_Dof_r - vBFxDof[j].Val[MAX_DIM*(k+1)+2] * Val_Dof_i ; Val[MAX_DIM*(k+1)+2] += vBFxDof[j].Val[MAX_DIM*k+2] * Val_Dof_i + vBFxDof[j].Val[MAX_DIM*(k+1)+2] * Val_Dof_r ; } } } *Type_Value = VECTOR ; break ; default : Msg(GERROR, "Unknown Form type in 'Pos_FemInterpolation'"); break; } } } if (Type_DefineQuantity != INTEGRALQUANTITY) break ; } /* while (1) ... */ if(UseNewGeo){ GeoData_P = (struct GeoData *)List_Pointer(GeoData_L, GeoDataNum); Geo_SetCurrentGeoData(Current.GeoData = GeoData_P) ; } GetDP_End ; }