is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.
| Cal_GalerkinTermOfFemEquation.cpp (getdp-3.4.0-source.tgz) | : | Cal_GalerkinTermOfFemEquation.cpp (getdp-3.5.0-source.tgz) |
|---|
|
| // GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege | | // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege |
| // | | // |
| // See the LICENSE.txt file for license information. Please report all | | // See the LICENSE.txt file for license information. Please report all |
| // issues on https://gitlab.onelab.info/getdp/getdp/issues. | | // issues on https://gitlab.onelab.info/getdp/getdp/issues. |
| // | | // |
| // Contributor(s): | | // Contributor(s): |
| // Johan Gyselinck | | // Johan Gyselinck |
| // Ruth Sabariego | | // Ruth Sabariego |
| // | | // |
| | |
| #include <map> | | #include <map> |
| | |
| skipping to change at line 30 | | skipping to change at line 30 |
|---|
| #include "Cal_AnalyticIntegration.h" | | #include "Cal_AnalyticIntegration.h" |
| #include "Cal_AssembleTerm.h" | | #include "Cal_AssembleTerm.h" |
| #include "Cal_GalerkinTermOfFemEquation.h" | | #include "Cal_GalerkinTermOfFemEquation.h" |
| #include "Get_DofOfElement.h" | | #include "Get_DofOfElement.h" |
| #include "Get_ElementSource.h" | | #include "Get_ElementSource.h" |
| #include "Get_Geometry.h" | | #include "Get_Geometry.h" |
| #include "Get_FunctionValue.h" | | #include "Get_FunctionValue.h" |
| #include "Pos_Search.h" | | #include "Pos_Search.h" |
| #include "Message.h" | | #include "Message.h" |
| | |
|
| extern struct Problem Problem_S ; | | extern struct Problem Problem_S; |
| extern struct CurrentData Current ; | | extern struct CurrentData Current; |
| | |
| std::map<int, bool> assDiag_done; | | std::map<int, bool> assDiag_done; |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l _ I n i t G a l e r k i n T e r m O f F e m E q u a t i o n */ | | /* C a l _ I n i t G a l e r k i n T e r m O f F e m E q u a t i o n */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Cal_InitGalerkinTermOfFemEquation(struct EquationTerm * EquationTerm_P, | | void Cal_InitGalerkinTermOfFemEquation( |
| struct QuantityStorage * QuantityStorage_ | | struct EquationTerm *EquationTerm_P, |
| P0, | | struct QuantityStorage *QuantityStorage_P0, |
| struct QuantityStorage * QuantityStorageN | | struct QuantityStorage *QuantityStorageNoDof, struct Dof *DofForNoDof_P) |
| oDof, | | |
| struct Dof * DofForNoDof_P) | | |
| { | | { |
|
| struct FemLocalTermActive * FI ; | | struct FemLocalTermActive *FI; |
| //extern int MH_Moving_Matrix_simple, MH_Moving_Matrix_probe, MH_Moving_Matrix | | extern int MHMoving_assemblyType; |
| _separate; | | |
| extern int MHMoving_assemblyType ; | | |
| | |
|
| FI = EquationTerm_P->Case.LocalTerm.Active ; | | FI = EquationTerm_P->Case.LocalTerm.Active; |
| | |
|
| FI->QuantityStorageEqu_P = QuantityStorage_P0 + | | FI->QuantityStorageEqu_P = |
| EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ; | | QuantityStorage_P0 + |
| | EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu; |
| | |
| Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu, | | Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu, |
|
| FI->QuantityStorageEqu_P, &FI->Type_FormEqu) ; | | FI->QuantityStorageEqu_P, &FI->Type_FormEqu); |
| | |
|
| if (EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof >= 0) { | | if(EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof >= 0) { |
| FI->QuantityStorageDof_P = QuantityStorage_P0 + | | FI->QuantityStorageDof_P = |
| EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof ; | | QuantityStorage_P0 + |
| FI->Type_DefineQuantityDof = FI->QuantityStorageDof_P->DefineQuantity->Type | | EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof; |
| ; | | FI->Type_DefineQuantityDof = FI->QuantityStorageDof_P->DefineQuantity->Type; |
| } | | } |
| else { | | else { |
|
| FI->QuantityStorageDof_P = QuantityStorageNoDof ; | | FI->QuantityStorageDof_P = QuantityStorageNoDof; |
| FI->Type_DefineQuantityDof = NODOF ; | | FI->Type_DefineQuantityDof = NODOF; |
| FI->DofForNoDof_P = DofForNoDof_P ; | | FI->DofForNoDof_P = DofForNoDof_P; |
| Dof_InitDofForNoDof(DofForNoDof_P, Current.NbrHar) ; | | Dof_InitDofForNoDof(DofForNoDof_P, Current.NbrHar); |
| QuantityStorageNoDof->BasisFunction[0].Dof = DofForNoDof_P ; | | QuantityStorageNoDof->BasisFunction[0].Dof = DofForNoDof_P; |
| } | | } |
| | |
| assDiag_done.clear(); | | assDiag_done.clear(); |
| | |
|
| ///*//kj+++ // brutal approach: comment this to avoid auto-symmetrization of Ja | | |
| cNL tensor with getdp (comment all this bloc) | | |
| // check for potentially symmetrical elementary matrix | | |
| FI->SymmetricalMatrix = | | FI->SymmetricalMatrix = |
| (EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu == | | (EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu == |
| EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof) && | | EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof) && |
| (EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu == | | (EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu == |
|
| EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof) ; | | EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof); |
| | |
|
| if(FI->SymmetricalMatrix){ | | if(FI->SymmetricalMatrix) { |
| // nonsymmetric if we play with test functions | | // nonsymmetric if we play with test functions |
|
| if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NON | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != |
| E){ | | CWQ_NONE) { |
| FI->SymmetricalMatrix = 0 ; | | FI->SymmetricalMatrix = 0; |
| } | | } |
| else{ | | else { |
| for(int i = 0; i < List_Nbr(EquationTerm_P->Case.LocalTerm.Term.WholeQuant | | for(int i = 0; |
| ity); i++){ | | i < List_Nbr(EquationTerm_P->Case.LocalTerm.Term.WholeQuantity); |
| struct WholeQuantity *WholeQuantity_P = (struct WholeQuantity*)List_Poin | | i++) { |
| ter | | struct WholeQuantity *WholeQuantity_P = |
| (EquationTerm_P->Case.LocalTerm.Term.WholeQuantity, i) ; | | (struct WholeQuantity *)List_Pointer( |
| | EquationTerm_P->Case.LocalTerm.Term.WholeQuantity, i); |
| // be on the safe side if we have a (noncommutative) vector product; | | // be on the safe side if we have a (noncommutative) vector product; |
|
| // FIXME: we should only do this id one of the arguments is a Dof | | // FIXME: we should only do this if one of the arguments is a Dof |
| if(WholeQuantity_P->Type == WQ_BINARYOPERATOR && | | if(WholeQuantity_P->Type == WQ_BINARYOPERATOR && |
| WholeQuantity_P->Case.Operator.TypeOperator == OP_CROSSPRODUCT) | | WholeQuantity_P->Case.Operator.TypeOperator == OP_CROSSPRODUCT) |
|
| FI->SymmetricalMatrix = 0 ; | | FI->SymmetricalMatrix = 0; |
| // FIXME: we should detect nonsymmetric tensors | | // FIXME: we should detect nonsymmetric tensors |
| } | | } |
| } | | } |
| } | | } |
|
| //*/ //kj+++ | | |
| //FI->SymmetricalMatrix = 0; //kj+++ // brutal approach: uncomment this to avoid | | |
| auto-symmetrization of JacNL tensor with getdp (uncomment this) | | |
| | |
|
| if (FI->SymmetricalMatrix) { | | if(FI->SymmetricalMatrix) { FI->Type_FormDof = FI->Type_FormEqu; } |
| FI->Type_FormDof = FI->Type_FormEqu ; | | |
| } | | |
| else { | | else { |
|
| switch (FI->Type_DefineQuantityDof) { | | switch(FI->Type_DefineQuantityDof) { |
| case LOCALQUANTITY : | | case LOCALQUANTITY: |
| Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof, | | Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof, |
|
| FI->QuantityStorageDof_P, &FI->Type_FormDof) ; | | FI->QuantityStorageDof_P, &FI->Type_FormDof); |
| break ; | | break; |
| case INTEGRALQUANTITY : | | case INTEGRALQUANTITY: |
| if(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof != NOOP){ | | if(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof != NOOP) { |
| Message::Error("No operator can act on an Integral Quantity"); | | Message::Error("No operator can act on an Integral Quantity"); |
| } | | } |
|
| FI->Type_FormDof = VECTOR ; /* we don't know the type a priori */ | | FI->Type_FormDof = VECTOR; /* we don't know the type a priori */ |
| FI->IntegralQuantityActive.IntegrationCase_L = | | FI->IntegralQuantityActive.IntegrationCase_L = |
|
| ((struct IntegrationMethod *) | | ((struct IntegrationMethod *)List_Pointer( |
| List_Pointer(Problem_S.IntegrationMethod, | | Problem_S.IntegrationMethod, |
| FI->QuantityStorageDof_P->DefineQuantity-> | | FI->QuantityStorageDof_P->DefineQuantity->IntegralQuantity |
| IntegralQuantity.IntegrationMethodIndex)) ->IntegrationCase | | .IntegrationMethodIndex)) |
| ; | | ->IntegrationCase; |
| FI->IntegralQuantityActive.CriterionIndex = | | FI->IntegralQuantityActive.CriterionIndex = |
|
| ((struct IntegrationMethod *) | | ((struct IntegrationMethod *)List_Pointer( |
| List_Pointer(Problem_S.IntegrationMethod, | | Problem_S.IntegrationMethod, |
| FI->QuantityStorageDof_P->DefineQuantity-> | | FI->QuantityStorageDof_P->DefineQuantity->IntegralQuantity |
| IntegralQuantity.IntegrationMethodIndex)) ->CriterionIndex | | .IntegrationMethodIndex)) |
| ; | | ->CriterionIndex; |
| FI->IntegralQuantityActive.JacobianCase_L = | | FI->IntegralQuantityActive.JacobianCase_L = |
|
| ((struct JacobianMethod *) | | ((struct JacobianMethod *)List_Pointer( |
| List_Pointer(Problem_S.JacobianMethod, | | Problem_S.JacobianMethod, FI->QuantityStorageDof_P->DefineQuantity |
| FI->QuantityStorageDof_P->DefineQuantity-> | | ->IntegralQuantity.JacobianMethodIndex)) |
| IntegralQuantity.JacobianMethodIndex)) ->JacobianCase ; | | ->JacobianCase; |
| break ; | | break; |
| case NODOF : | | case NODOF: FI->Type_FormDof = SCALAR; break; |
| FI->Type_FormDof = SCALAR ; | | |
| break ; | | |
| } | | } |
| } | | } |
| | |
| FI->Type_ValueDof = Get_ValueFromForm(FI->Type_FormDof); | | FI->Type_ValueDof = Get_ValueFromForm(FI->Type_FormDof); |
| | |
| /* G e t I n t e g r a t i o n M e t h o d */ | | /* G e t I n t e g r a t i o n M e t h o d */ |
| /* -------------------------------------------- */ | | /* -------------------------------------------- */ |
| | |
|
| if(EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex < 0){ | | if(EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex < 0) { |
| Message::Error("Integration method missing in equation term"); | | Message::Error("Integration method missing in equation term"); |
| FI->IntegrationCase_L = 0; | | FI->IntegrationCase_L = 0; |
| } | | } |
|
| else{ | | else { |
| FI->IntegrationCase_L = | | FI->IntegrationCase_L = |
|
| ((struct IntegrationMethod *) | | ((struct IntegrationMethod *)List_Pointer( |
| List_Pointer(Problem_S.IntegrationMethod, | | Problem_S.IntegrationMethod, |
| EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex)) | | EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex)) |
| ->IntegrationCase ; | | ->IntegrationCase; |
| | |
| FI->CriterionIndex = | | FI->CriterionIndex = |
|
| ((struct IntegrationMethod *) | | ((struct IntegrationMethod *)List_Pointer( |
| List_Pointer(Problem_S.IntegrationMethod, | | Problem_S.IntegrationMethod, |
| EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex)) | | EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex)) |
| ->CriterionIndex ; | | ->CriterionIndex; |
| } | | } |
| | |
| /* G e t J a c o b i a n M e t h o d */ | | /* G e t J a c o b i a n M e t h o d */ |
| /* -------------------------------------- */ | | /* -------------------------------------- */ |
| | |
|
| if(EquationTerm_P->Case.LocalTerm.JacobianMethodIndex < 0){ | | if(EquationTerm_P->Case.LocalTerm.JacobianMethodIndex < 0) { |
| Message::Error("Jacobian method missing in equation term"); | | Message::Error("Jacobian method missing in equation term"); |
| FI->JacobianCase_L = 0; | | FI->JacobianCase_L = 0; |
| } | | } |
|
| else{ | | else { |
| FI->JacobianCase_L = | | FI->JacobianCase_L = ((struct JacobianMethod *)List_Pointer( |
| ((struct JacobianMethod *) | | Problem_S.JacobianMethod, |
| List_Pointer(Problem_S.JacobianMethod, | | EquationTerm_P->Case.LocalTerm.JacobianMethodIndex)) |
| EquationTerm_P->Case.LocalTerm.JacobianMethodIndex)) | | ->JacobianCase; |
| ->JacobianCase ; | | |
| | |
| FI->JacobianCase_P0 = | | FI->JacobianCase_P0 = |
| (FI->NbrJacobianCase = List_Nbr(FI->JacobianCase_L)) ? | | (FI->NbrJacobianCase = List_Nbr(FI->JacobianCase_L)) ? |
|
| (struct JacobianCase*)List_Pointer(FI->JacobianCase_L, 0) : NULL ; | | (struct JacobianCase *)List_Pointer(FI->JacobianCase_L, 0) : |
| | NULL; |
| } | | } |
| | |
| FI->Flag_ChangeCoord = | | FI->Flag_ChangeCoord = |
|
| ( FI->SymmetricalMatrix || | | (FI->SymmetricalMatrix || |
| !( | | !(((FI->Type_FormEqu == FORM0 || FI->Type_FormEqu == FORM0P) && |
| ( (FI->Type_FormEqu == FORM0 || FI->Type_FormEqu == FORM0P) && | | (FI->Type_FormDof == FORM3 || FI->Type_FormDof == FORM3P)) || |
| (FI->Type_FormDof == FORM3 || FI->Type_FormDof == FORM3P) ) || | | /* |
| /* | | ( (FI->Type_FormEqu == FORM1 || FI->Type_FormEqu == FORM1P) && |
| ( (FI->Type_FormEqu == FORM1 || FI->Type_FormEqu == FORM1P) && | | (FI->Type_FormDof == FORM2 || FI->Type_FormDof == FORM2P) ) || |
| (FI->Type_FormDof == FORM2 || FI->Type_FormDof == FORM2P) ) || | | ( (FI->Type_FormEqu == FORM2 || FI->Type_FormEqu == FORM2P) && |
| ( (FI->Type_FormEqu == FORM2 || FI->Type_FormEqu == FORM2P) && | | (FI->Type_FormDof == FORM1 || FI->Type_FormDof == FORM1P) ) || |
| (FI->Type_FormDof == FORM1 || FI->Type_FormDof == FORM1P) ) || | | */ |
| */ | | ((FI->Type_FormEqu == FORM3 || FI->Type_FormEqu == FORM3P) && |
| ( (FI->Type_FormEqu == FORM3 || FI->Type_FormEqu == FORM3P) && | | (FI->Type_FormDof == FORM0 || |
| (FI->Type_FormDof == FORM0 || FI->Type_FormDof == FORM0P) ) | | FI->Type_FormDof == |
| ) | | FORM0P)))) || /* For operators on VECTOR's (To be extended) */ |
| ) | | (FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR) || |
| || /* For operators on VECTOR's (To be extended) */ | | (FI->Type_DefineQuantityDof == INTEGRALQUANTITY); |
| (FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR) | | |
| || | | if(FI->Flag_ChangeCoord) { |
| (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) ; | | FI->Flag_InvJac = |
| | ((FI->Type_FormEqu == FORM1) || |
| if (FI->Flag_ChangeCoord){ | | (!FI->SymmetricalMatrix && (FI->Type_FormDof == FORM1)) || |
| FI->Flag_InvJac = ( (FI->Type_FormEqu == FORM1) || | | /* For operators on VECTOR's (To be extended) */ |
| (!FI->SymmetricalMatrix && (FI->Type_FormDof == FORM1)) | | (FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR) || |
| || | | (EquationTerm_P->Case.LocalTerm.Term.QuantityIndexPost == 1)); |
| /* For operators on VECTOR's (To be extended) */ | | |
| (FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTO | | |
| R) || | | |
| (EquationTerm_P->Case.LocalTerm.Term.QuantityIndexPost = | | |
| = 1) ) ; | | |
| } | | } |
| | |
| /* G e t C h a n g e O f C o o r d i n a t e s */ | | /* G e t C h a n g e O f C o o r d i n a t e s */ |
| /* ---------------------------------------------- */ | | /* ---------------------------------------------- */ |
| | |
| FI->xChangeOfCoordinatesEqu = | | FI->xChangeOfCoordinatesEqu = |
|
| (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormEqu)
; | | (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormEqu); |
| | |
|
| if (!FI->SymmetricalMatrix) | | if(!FI->SymmetricalMatrix) |
| FI->xChangeOfCoordinatesDof = | | FI->xChangeOfCoordinatesDof = (void (*)())Get_ChangeOfCoordinates( |
| (void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormDof | | FI->Flag_ChangeCoord, FI->Type_FormDof); |
| ) ; | | |
| else | | else |
|
| FI->xChangeOfCoordinatesDof = | | FI->xChangeOfCoordinatesDof = (void (*)())Get_ChangeOfCoordinates( |
| (void (*)())Get_ChangeOfCoordinates(0, FI->Type_FormDof) ; /* Used only fo | | 0, FI->Type_FormDof); /* Used only for transfer */ |
| r transfer */ | | |
| | |
| /* G e t C a l _ P r o d u c t x */ | | /* G e t C a l _ P r o d u c t x */ |
| /* -------------------------------- */ | | /* -------------------------------- */ |
| | |
|
| switch (FI->Type_FormEqu) { | | switch(FI->Type_FormEqu) { |
| case FORM1 : case FORM1S : | | case FORM1: |
| case FORM2 : case FORM2P : case FORM2S : | | case FORM1S: |
| case VECTOR : | | case FORM2: |
| FI->Cal_Productx = (double (*)())Cal_Product123 ; break ; | | case FORM2P: |
| case FORM1P : | | case FORM2S: |
| case VECTORP : | | case VECTOR: FI->Cal_Productx = (double (*)())Cal_Product123; break; |
| FI->Cal_Productx = (double (*)())Cal_Product3 ; break ; | | case FORM1P: |
| case FORM0 : | | case VECTORP: FI->Cal_Productx = (double (*)())Cal_Product3; break; |
| case FORM3 : case FORM3P : | | case FORM0: |
| case SCALAR : | | case FORM3: |
| FI->Cal_Productx = (double (*)())Cal_Product1 ; break ; | | case FORM3P: |
| default : | | case SCALAR: FI->Cal_Productx = (double (*)())Cal_Product1; break; |
| | default: |
| Message::Error("Unknown type of Form (%d)", FI->Type_FormEqu); | | Message::Error("Unknown type of Form (%d)", FI->Type_FormEqu); |
|
| FI->Cal_Productx = (double (*)())Cal_Product123 ; break ; | | FI->Cal_Productx = (double (*)())Cal_Product123; |
| | break; |
| } | | } |
| | |
| /* G e t F u n c t i o n _ A s s e m b l e T e r m */ | | /* G e t F u n c t i o n _ A s s e m b l e T e r m */ |
| /* ------------------------------------------------- */ | | /* ------------------------------------------------- */ |
| | |
|
| switch (EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative) { | | switch(EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative) { |
| case NODT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | case NODT_: |
| _NoDt ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NoDt; |
| case DT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | break; |
| _Dt ; break; | | case DT_: FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_Dt; break; |
| case DTDOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | case DTDOF_: |
| _DtDof ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDof; |
| case DTDT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | break; |
| _DtDt ; break; | | case DTDT_: |
| case DTDTDOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDt; |
| _DtDtDof ; break; | | break; |
| case DTDTDTDOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | case DTDTDOF_: |
| _DtDtDtDof ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDtDof; |
| case DTDTDTDTDOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | break; |
| _DtDtDtDtDof ; break; | | case DTDTDTDOF_: |
| case DTDTDTDTDTDOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDtDtDof; |
| _DtDtDtDtDtDof; break; | | break; |
| case JACNL_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | case DTDTDTDTDOF_: |
| _JacNL ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDtDtDtDof; |
| case DTDOFJACNL_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | break; |
| _DtDofJacNL ; break; | | case DTDTDTDTDTDOF_: |
| case NEVERDT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDtDtDtDtDof; |
| _NeverDt ; break; | | break; |
| case DTNL_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm | | case JACNL_: |
| _DtNL ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_JacNL; |
| | break; |
| | case DTDOFJACNL_: |
| | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDofJacNL; |
| | break; |
| | case NEVERDT_: |
| | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NeverDt; |
| | break; |
| | case DTNL_: |
| | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtNL; |
| | break; |
| // nleigchange | | // nleigchange |
|
| case NLEIG1DOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTer | | case NLEIG1DOF_: |
| m_NLEig1Dof ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NLEig1Dof; |
| case NLEIG2DOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTer | | break; |
| m_NLEig2Dof ; break; | | case NLEIG2DOF_: |
| case NLEIG3DOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTer | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NLEig2Dof; |
| m_NLEig3Dof ; break; | | break; |
| case NLEIG4DOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTer | | case NLEIG3DOF_: |
| m_NLEig4Dof ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NLEig3Dof; |
| case NLEIG5DOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTer | | break; |
| m_NLEig5Dof ; break; | | case NLEIG4DOF_: |
| case NLEIG6DOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTer | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NLEig4Dof; |
| m_NLEig6Dof ; break; | | break; |
| | case NLEIG5DOF_: |
| | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NLEig5Dof; |
| | break; |
| | case NLEIG6DOF_: |
| | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NLEig6Dof; |
| | break; |
| | |
|
| default : | | default: |
| Message::Error("Unknown type of Operator for Galerkin term (%d)", | | Message::Error("Unknown type of Operator for Galerkin term (%d)", |
| EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative); | | EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative); |
|
| FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NoDt ; break; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NoDt; |
| | break; |
| } | | } |
| | |
| if(MHMoving_assemblyType) | | if(MHMoving_assemblyType) |
| FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MHMoving; | | FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MHMoving; |
| | |
|
| /* | | |
| if (MH_Moving_Matrix_simple) { | | |
| FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_simple ; | | |
| } | | |
| if (MH_Moving_Matrix_probe) { | | |
| FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_probe ; | | |
| } | | |
| if (MH_Moving_Matrix_separate) { | | |
| FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_separate | | |
| ; | | |
| } | | |
| */ | | |
| | | |
| // TODO: if JACNL_, say to Cal_Init to assemble later in Jac, otherwise | | // TODO: if JACNL_, say to Cal_Init to assemble later in Jac, otherwise |
| // assemble in the system matrix | | // assemble in the system matrix |
| | |
|
| /* initialisation of MHBilinear-term (nonlinear multi-harmonics) if necessary
*/ | | // initialisation of MHBilinear-term (nonlinear multi-harmonics) |
| Cal_InitGalerkinTermOfFemEquation_MHBilinear(EquationTerm_P); | | Cal_InitGalerkinTermOfFemEquation_MHBilinear(EquationTerm_P); |
| | |
|
| /* Full_Matrix */ | | if(EquationTerm_P->Case.LocalTerm.Full_Matrix) { |
| if (EquationTerm_P->Case.LocalTerm.Full_Matrix) { | | |
| FI->Full_Matrix = 1; | | FI->Full_Matrix = 1; |
|
| FI->FirstElements = List_Create(20, 10, sizeof(struct FirstElement)) ; | | FI->FirstElements = List_Create(20, 10, sizeof(struct FirstElement)); |
| } | | } |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l _ E n d G a l e r k i n T e r m O f F e m E q u a t i o n */ | | /* C a l _ E n d G a l e r k i n T e r m O f F e m E q u a t i o n */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Cal_EndGalerkinTermOfFemEquation() | | void Cal_EndGalerkinTermOfFemEquation() { assDiag_done.clear(); } |
| { | | |
| assDiag_done.clear(); | | |
| } | | |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l _ a p p l y M e t r i c T e n s o r */ | | /* C a l _ a p p l y M e t r i c T e n s o r */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Cal_applyMetricTensor(struct EquationTerm * EquationTerm_P, | | void Cal_applyMetricTensor(struct EquationTerm *EquationTerm_P, |
| struct FemLocalTermActive * FI, | | struct FemLocalTermActive *FI, |
| struct QuantityStorage * QuantityStorage_P0, | | struct QuantityStorage *QuantityStorage_P0, |
| int Nbr_Dof, | | int Nbr_Dof, struct Value vBFxDof[]) |
| struct Value vBFxDof[]) | | |
| { | | { |
|
| int j; | | int j; |
| int mi; | | int mi; |
| struct Value S; | | struct Value S; |
| struct Value detS; | | struct Value detS; |
| | |
| mi = EquationTerm_P->Case.LocalTerm.ExpressionIndexForMetricTensor; | | mi = EquationTerm_P->Case.LocalTerm.ExpressionIndexForMetricTensor; |
| if(mi == -1) return; | | if(mi == -1) return; |
|
| Get_ValueOfExpression | | Get_ValueOfExpression( |
| ((struct Expression*)List_Pointer(Problem_S.Expression, mi), | | (struct Expression *)List_Pointer(Problem_S.Expression, mi), |
| QuantityStorage_P0, Current.u, Current.v, Current.w, &S) ; | | QuantityStorage_P0, Current.u, Current.v, Current.w, &S); |
| | |
| if(S.Type == SCALAR) { | | if(S.Type == SCALAR) { |
| S.Type = TENSOR_DIAG; | | S.Type = TENSOR_DIAG; |
| S.Val[1] = S.Val[0]; | | S.Val[1] = S.Val[0]; |
| S.Val[2] = S.Val[0]; | | S.Val[2] = S.Val[0]; |
| } | | } |
| if(S.Type != TENSOR_SYM && S.Type != TENSOR && S.Type != TENSOR_DIAG) { | | if(S.Type != TENSOR_SYM && S.Type != TENSOR && S.Type != TENSOR_DIAG) { |
| Message::Error("Cannot interpret field type %s as metric tensor", | | Message::Error("Cannot interpret field type %s as metric tensor", |
| Get_StringForDefine(Field_Type, S.Type)); | | Get_StringForDefine(Field_Type, S.Type)); |
| return; | | return; |
| } | | } |
| | |
| Cal_DetValue(&S, &detS); | | Cal_DetValue(&S, &detS); |
| detS.Val[0] = sqrt(fabs(detS.Val[0])); | | detS.Val[0] = sqrt(fabs(detS.Val[0])); |
| | |
|
| switch (FI->Type_FormDof) { | | switch(FI->Type_FormDof) { |
| case FORM1 : case FORM1S : case FORM1P : | | case FORM1: |
| | case FORM1S: |
| | case FORM1P: |
| Cal_InvertValue(&S, &S); | | Cal_InvertValue(&S, &S); |
|
| for (j = 0 ; j < Nbr_Dof ; j++) { | | for(j = 0; j < Nbr_Dof; j++) { |
| Cal_ProductValue(&S, &vBFxDof[j], &vBFxDof[j]); | | Cal_ProductValue(&S, &vBFxDof[j], &vBFxDof[j]); |
| Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); | | Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); |
| } | | } |
| break; | | break; |
|
| case FORM2 : case FORM2S : case FORM2P : | | case FORM2: |
| | case FORM2S: |
| | case FORM2P: |
| Cal_InvertValue(&detS, &detS); | | Cal_InvertValue(&detS, &detS); |
|
| for (j = 0 ; j < Nbr_Dof ; j++) { | | for(j = 0; j < Nbr_Dof; j++) { |
| Cal_ProductValue(&S, &vBFxDof[j], &vBFxDof[j]); | | Cal_ProductValue(&S, &vBFxDof[j], &vBFxDof[j]); |
| Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); | | Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); |
| } | | } |
| break; | | break; |
|
| case FORM3 : case FORM3S : case FORM3P : | | case FORM3: |
| | case FORM3S: |
| | case FORM3P: |
| Cal_InvertValue(&detS, &detS); | | Cal_InvertValue(&detS, &detS); |
|
| for (j = 0 ; j < Nbr_Dof ; j++) { | | for(j = 0; j < Nbr_Dof; j++) { |
| Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); | | Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); |
| } | | } |
| break; | | break; |
|
| case FORM0 : case FORM0S : case FORM0P : | | case FORM0: |
| for (j = 0 ; j < Nbr_Dof ; j++) { | | case FORM0S: |
| | case FORM0P: |
| | for(j = 0; j < Nbr_Dof; j++) { |
| Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); | | Cal_ProductValue(&detS, &vBFxDof[j], &vBFxDof[j]); |
| } | | } |
| break; | | break; |
| default: | | default: |
| Message::Error("Cannot use metric tensor with field type %s", | | Message::Error("Cannot use metric tensor with field type %s", |
| Get_StringForDefine(Field_Type, FI->Type_FormDof)); | | Get_StringForDefine(Field_Type, FI->Type_FormDof)); |
| break; | | break; |
| } | | } |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l _ v B F x D o f */ | | /* C a l _ v B F x D o f */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Cal_vBFxDof(struct EquationTerm * EquationTerm_P, | | int Cal_vBFxDof(struct EquationTerm *EquationTerm_P, |
| struct FemLocalTermActive * FI, | | struct FemLocalTermActive *FI, |
| struct QuantityStorage * QuantityStorage_P0, | | struct QuantityStorage *QuantityStorage_P0, |
| struct QuantityStorage * QuantityStorageDof_P, | | struct QuantityStorage *QuantityStorageDof_P, int Nbr_Dof, |
| int Nbr_Dof, | | void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])( |
| void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS]) | | struct Element *Element, int NumEntity, double u, double v, |
| (struct Element * Element, int NumEntity, | | double w, double Value[]), |
| double u, double v, double w, double Value[]), | | double vBFxEqu[][MAX_DIM], struct Value vBFxDof[]) |
| double vBFxEqu[][MAX_DIM], | | |
| struct Value vBFxDof[]) | | |
| { | | { |
|
| double vBFuDof[NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ; | | double vBFuDof[NBR_MAX_BASISFUNCTIONS][MAX_DIM]; |
| double u, v, w ; | | double u, v, w; |
| struct Value CoefPhys ; | | struct Value CoefPhys; |
| struct Element *E ; | | struct Element *E; |
| int i, j ; | | int i, j; |
| | |
| if(EquationTerm_P->Case.LocalTerm.Term.DofInTrace){ | | |
| E = Current.Element->ElementTrace ; | | |
| Current.x = Current.y = Current.z = 0. ; | | |
| for (i = 0 ; i < Current.Element->GeoElement->NbrNodes ; i++) { | | |
| Current.x += Current.Element->x[i] * Current.Element->n[i] ; | | |
| Current.y += Current.Element->y[i] * Current.Element->n[i] ; | | |
| Current.z += Current.Element->z[i] * Current.Element->n[i] ; | | |
| } | | |
| xyz2uvwInAnElement(E, Current.x, Current.y, Current.z, | | |
| &Current.ut, &Current.vt, &Current.wt) ; | | |
| u = Current.ut ; | | |
| v = Current.vt ; | | |
| w = Current.wt ; | | |
| } | | |
| else{ | | |
| E = Current.Element ; | | |
| u = Current.u ; | | |
| v = Current.v ; | | |
| w = Current.w ; | | |
| } | | |
| | |
| // initialize vBFxDof to zero; this allows to perform e.g. [0, {d u}] without | | // initialize vBFxDof to zero; this allows to perform e.g. [0, {d u}] without |
| // having to explicitly use [Vector[0,0,0], {d u}] ; if this is too slow, we | | // having to explicitly use [Vector[0,0,0], {d u}] ; if this is too slow, we |
| // should check vBFxDof[j].Type against FI->Type_FormEqu before calling | | // should check vBFxDof[j].Type against FI->Type_FormEqu before calling |
| // FI->Cal_Productx to report errors | | // FI->Cal_Productx to report errors |
|
| for (j = 0 ; j < Nbr_Dof ; j++) | | for(j = 0; j < Nbr_Dof; j++) Cal_ZeroValue(&vBFxDof[j]); |
| Cal_ZeroValue(&vBFxDof[j]); | | |
| | |
|
| // shape functions, integral quantity or dummy | | if(EquationTerm_P->Case.LocalTerm.Term.DofInTrace) { |
| | E = Current.Element->ElementTrace; |
| | Current.x = Current.y = Current.z = 0.; |
| | for(i = 0; i < Current.Element->GeoElement->NbrNodes; i++) { |
| | Current.x += Current.Element->x[i] * Current.Element->n[i]; |
| | Current.y += Current.Element->y[i] * Current.Element->n[i]; |
| | Current.z += Current.Element->z[i] * Current.Element->n[i]; |
| | } |
| | // TODO we might want to make this a parameter |
| | double tol = Current.GeoData->CharacteristicLength * 1.e-4; |
| | if(!PointInElement(E, nullptr, Current.x, Current.y, Current.z, &Current.ut, |
| | &Current.vt, &Current.wt, tol)) { |
| | return 0; // we're done, no contribution from this Trace element |
| | } |
| | u = Current.ut; |
| | v = Current.vt; |
| | w = Current.wt; |
| | } |
| | else { |
| | E = Current.Element; |
| | u = Current.u; |
| | v = Current.v; |
| | w = Current.w; |
| | } |
| | |
|
| if (!FI->SymmetricalMatrix) { | | // shape functions, integral quantity or dummy |
| | |
|
| switch (FI->Type_DefineQuantityDof) { | | if(!FI->SymmetricalMatrix) { |
| case LOCALQUANTITY : | | switch(FI->Type_DefineQuantityDof) { |
| for (j = 0 ; j < Nbr_Dof ; j++) { | | case LOCALQUANTITY: |
| xFunctionBFDof[j] | | for(j = 0; j < Nbr_Dof; j++) { |
| (E, | | xFunctionBFDof[j]( |
| QuantityStorageDof_P->BasisFunction[j].NumEntityInElement+1, | | E, QuantityStorageDof_P->BasisFunction[j].NumEntityInElement + 1, u, |
| u, v, w, vBFuDof[j]) ; | | v, w, vBFuDof[j]); |
| ((void (*)(struct Element*, double*, double*)) | | ((void (*)(struct Element *, double *, |
| FI->xChangeOfCoordinatesDof) (E, vBFuDof[j], vBFxDof[j].Val) ; | | double *))FI->xChangeOfCoordinatesDof)(E, vBFuDof[j], |
| vBFxDof[j].Type = FI->Type_ValueDof ; | | vBFxDof[j].Val); |
| if(Current.NbrHar > 1) Cal_SetHarmonicValue(&vBFxDof[j]) ; | | vBFxDof[j].Type = FI->Type_ValueDof; |
| | | if(Current.NbrHar > 1) Cal_SetHarmonicValue(&vBFxDof[j]); |
| /* temp (rather add QuantityStorage_P to CurrentData) */ | | |
| Current.NumEntities[j] = | | /* temp (rather add QuantityStorage_P to CurrentData) */ |
| QuantityStorageDof_P->BasisFunction[j].CodeEntity; | | Current.NumEntities[j] = |
| } | | QuantityStorageDof_P->BasisFunction[j].CodeEntity; |
| break ; | | } |
| case INTEGRALQUANTITY : | | break; |
| if (FI->IntegralQuantityActive.IntegrationCase_P->Type == ANALYTIC) | | case INTEGRALQUANTITY: |
| Cal_AnalyticIntegralQuantity (Current.Element, | | if(FI->IntegralQuantityActive.IntegrationCase_P->Type == ANALYTIC) |
| QuantityStorageDof_P, Nbr_Dof, | | Cal_AnalyticIntegralQuantity(Current.Element, QuantityStorageDof_P, |
| (void (**)())xFunctionBFDof, vBFxDof) ; | | Nbr_Dof, (void (**)())xFunctionBFDof, |
| | vBFxDof); |
| | |
| else | | else |
|
| Cal_NumericalIntegralQuantity (Current.Element, | | Cal_NumericalIntegralQuantity( |
| &FI->IntegralQuantityActive, | | Current.Element, &FI->IntegralQuantityActive, QuantityStorage_P0, |
| QuantityStorage_P0, QuantityStorageDof_P, | | QuantityStorageDof_P, FI->Type_DefineQuantityDof, Nbr_Dof, |
| FI->Type_DefineQuantityDof, Nbr_Dof, | | (void (**)())xFunctionBFDof, vBFxDof); |
| (void (**)())xFunctionBFDof, vBFxDof) ; | | FI->Type_ValueDof = FI->Type_FormDof = |
| FI->Type_ValueDof = FI->Type_FormDof = vBFxDof[0].Type; /* now this type i | | vBFxDof[0].Type; /* now this type is correct */ |
| s correct */ | | break; |
| break ; | | case NODOF: /* Supprimer le DofForNoDof_P de la structure dans Data_Active.h |
| case NODOF : /* Supprimer le DofForNoDof_P de la structure dans Data_Active | | */ |
| .h */ | | /* QuantityStorageDof_P->BasisFunction[0].Dof = FI->DofForNoDof_P ; |
| /* QuantityStorageDof_P->BasisFunction[0].Dof = FI->DofForNoDof_P ; * | | */ |
| / | | break; |
| break ; | | |
| } | | } |
| } | | } |
| else { | | else { |
|
| for (j = 0 ; j < Nbr_Dof ; j++){ | | for(j = 0; j < Nbr_Dof; j++) { |
| ((void (*)(struct Element*, double*, double*)) | | ((void (*)(struct Element *, double *, |
| FI->xChangeOfCoordinatesDof) (Current.Element, vBFxEqu[j], vBFxDof[j].Val | | double *))FI->xChangeOfCoordinatesDof)( |
| ) ; | | Current.Element, vBFxEqu[j], vBFxDof[j].Val); |
| vBFxDof[j].Type = FI->Type_ValueDof ; | | vBFxDof[j].Type = FI->Type_ValueDof; |
| if(Current.NbrHar > 1) Cal_SetHarmonicValue(&vBFxDof[j]) ; | | if(Current.NbrHar > 1) Cal_SetHarmonicValue(&vBFxDof[j]); |
| } | | } |
| } | | } |
| | |
| /* Compute remaining factors in the term */ | | /* Compute remaining factors in the term */ |
| | |
|
| if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == CWQ_DOF) {} |
| CWQ_DOF) { | | else if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == |
| } | | CWQ_EXP_TIME_DOF) { |
| else if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == | | Get_ValueOfExpression( |
| CWQ_EXP_TIME_DOF) { | | (struct Expression *)List_Pointer( |
| Get_ValueOfExpression | | Problem_S.Expression, |
| ((struct Expression*)List_Pointer | | EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical), |
| (Problem_S.Expression, | | QuantityStorage_P0, Current.u, Current.v, Current.w, &CoefPhys); |
| EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical), | | for(j = 0; j < Nbr_Dof; j++) |
| QuantityStorage_P0, Current.u, Current.v, Current.w, | | Cal_ProductValue(&CoefPhys, &vBFxDof[j], &vBFxDof[j]); |
| &CoefPhys) ; | | |
| for (j = 0 ; j < Nbr_Dof ; j++) | | |
| Cal_ProductValue(&CoefPhys, &vBFxDof[j], &vBFxDof[j]) ; | | |
| } | | } |
| else | | else |
|
| Cal_WholeQuantity | | Cal_WholeQuantity( |
| (Current.Element, QuantityStorage_P0, | | Current.Element, QuantityStorage_P0, |
| EquationTerm_P->Case.LocalTerm.Term.WholeQuantity, | | EquationTerm_P->Case.LocalTerm.Term.WholeQuantity, Current.u, Current.v, |
| Current.u, Current.v, Current.w, | | Current.w, EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity, |
| EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity, | | Nbr_Dof, vBFxDof); |
| Nbr_Dof, vBFxDof) ; | | |
| | |
| /* Compute using metric tensor if defined */ | | /* Compute using metric tensor if defined */ |
|
| Cal_applyMetricTensor(EquationTerm_P, FI, QuantityStorage_P0, | | Cal_applyMetricTensor(EquationTerm_P, FI, QuantityStorage_P0, Nbr_Dof, |
| Nbr_Dof, vBFxDof); | | vBFxDof); |
| | |
|
| | return 1; |
| } | | } |
| | |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| /* C a l _ T e r m O f F e m E q u a t i o n */ | | /* C a l _ T e r m O f F e m E q u a t i o n */ |
| /* ------------------------------------------------------------------------ */ | | /* ------------------------------------------------------------------------ */ |
| | |
|
| void Cal_GalerkinTermOfFemEquation(struct Element * Element, | | void Cal_GalerkinTermOfFemEquation(struct Element *Element, |
| struct EquationTerm * EquationTerm_P, | | struct EquationTerm *EquationTerm_P, |
| struct QuantityStorage * QuantityStorage_P0) | | struct QuantityStorage *QuantityStorage_P0) |
| { | | { |
|
| struct FemLocalTermActive * FI ; | | struct FemLocalTermActive *FI; |
| struct QuantityStorage * QuantityStorageEqu_P, * QuantityStorageDof_P ; | | struct QuantityStorage *QuantityStorageEqu_P, *QuantityStorageDof_P; |
| struct IntegrationCase * IntegrationCase_P ; | | struct IntegrationCase *IntegrationCase_P; |
| struct Quadrature * Quadrature_P ; | | struct Quadrature *Quadrature_P; |
| struct Value vBFxDof [NBR_MAX_BASISFUNCTIONS], CoefPhys ; | | struct Value vBFxDof[NBR_MAX_BASISFUNCTIONS], CoefPhys; |
| struct Value CanonicExpression_Equ, V1, V2; | | struct Value CanonicExpression_Equ, V1, V2; |
| | | |
| int Nbr_Equ, Nbr_Dof = 0; | | int Nbr_Equ, Nbr_Dof = 0; |
| int i, j, k, Type_Dimension, Nbr_IntPoints, i_IntPoint ; | | int i, j, k, Type_Dimension, Nbr_IntPoints, i_IntPoint; |
| int NextElement ; | | int NextElement; |
| | | |
| double weight, Factor = 1. ; | | double weight, Factor = 1.; |
| double vBFuEqu [NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ; | | double vBFuEqu[NBR_MAX_BASISFUNCTIONS][MAX_DIM]; |
| double vBFxEqu [NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ; | | double vBFxEqu[NBR_MAX_BASISFUNCTIONS][MAX_DIM]; |
| double Ek [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_HARMONIC | | double Ek[NBR_MAX_BASISFUNCTIONS][NBR_MAX_BASISFUNCTIONS][NBR_MAX_HARMONIC]; |
| ] ; | | |
| | | void (*xFunctionBFEqu[NBR_MAX_BASISFUNCTIONS])( |
| void (*xFunctionBFEqu[NBR_MAX_BASISFUNCTIONS]) | | struct Element * Element, int NumEntity, double u, double v, double w, |
| (struct Element * Element, int NumEntity, | | double Value[]); |
| double u, double v, double w, double Value[] ) ; | | void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])( |
| void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS]) | | struct Element * Element, int NumEntity, double u, double v, double w, |
| (struct Element * Element, int NumEntity, | | double Value[]); |
| double u, double v, double w, double Value[] ) ; | | double (*Get_Jacobian)(struct Element *, MATRIX3x3 *); |
| double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ; | | void (*Get_IntPoint)(int, int, double *, double *, double *, double *); |
| void (*Get_IntPoint)(int,int,double*,double*,double*,double*); | | |
| | |
| extern int Flag_RHS; | | extern int Flag_RHS; |
| | |
| Current.flagAssDiag = 0; /*+++prov*/ | | Current.flagAssDiag = 0; /*+++prov*/ |
| | |
|
| FI = EquationTerm_P->Case.LocalTerm.Active ; | | FI = EquationTerm_P->Case.LocalTerm.Active; |
| | |
| /* treatment of MHBilinear-term in separate routine */ | | /* treatment of MHBilinear-term in separate routine */ |
|
| if (FI->MHBilinear) { | | if(FI->MHBilinear) { |
| /* if only the RHS of the system is to be calculated | | /* if only the RHS of the system is to be calculated |
| (in case of adaptive relaxation of the Newton-Raphson scheme) | | (in case of adaptive relaxation of the Newton-Raphson scheme) |
| the (expensive and redundant) calculation of this term can be skipped */ | | the (expensive and redundant) calculation of this term can be skipped */ |
|
| if (!Flag_RHS) | | if(!Flag_RHS) |
| Cal_GalerkinTermOfFemEquation_MHBilinear(Element, EquationTerm_P, Quantity | | Cal_GalerkinTermOfFemEquation_MHBilinear(Element, EquationTerm_P, |
| Storage_P0) ; | | QuantityStorage_P0); |
| return; | | return; |
| } | | } |
| | |
|
| QuantityStorageEqu_P = FI->QuantityStorageEqu_P ; | | QuantityStorageEqu_P = FI->QuantityStorageEqu_P; |
| QuantityStorageDof_P = FI->QuantityStorageDof_P ; | | QuantityStorageDof_P = FI->QuantityStorageDof_P; |
| | |
| /* skip computation completely if term does not contribute to RHS. This is OK, | | /* skip computation completely if term does not contribute to RHS. This is OK, |
| but the speed-up is not the best, due to the time-consuming--but | | but the speed-up is not the best, due to the time-consuming--but |
| necessary-- initializations that still need to be done before arriving at | | necessary-- initializations that still need to be done before arriving at |
| this point in the assembly process. For best performance use | | this point in the assembly process. For best performance use |
| GenerateRHSGroup instead of GenerateRHS, and include any RHS-contributions | | GenerateRHSGroup instead of GenerateRHS, and include any RHS-contributions |
| (elements containing fixed dofs or non-dof expressions) in the given | | (elements containing fixed dofs or non-dof expressions) in the given |
| groups */ | | groups */ |
|
| if(Current.DofData->Flag_RHS){ | | if(Current.DofData->Flag_RHS) { |
| if(FI->Type_DefineQuantityDof == LOCALQUANTITY){ | | if(FI->Type_DefineQuantityDof == LOCALQUANTITY) { |
| bool skip = true; | | bool skip = true; |
|
| int Nbr_Dof = FI->SymmetricalMatrix ? QuantityStorageEqu_P->NbrElementaryB | | int Nbr_Dof = FI->SymmetricalMatrix ? |
| asisFunction : | | QuantityStorageEqu_P->NbrElementaryBasisFunction : |
| QuantityStorageDof_P->NbrElementaryBasisFunction; | | QuantityStorageDof_P->NbrElementaryBasisFunction; |
| for (int j = 0 ; j < Nbr_Dof ; j++){ | | for(int j = 0; j < Nbr_Dof; j++) { |
| Dof *Dof_P = QuantityStorageDof_P->BasisFunction[j].Dof; | | Dof *Dof_P = QuantityStorageDof_P->BasisFunction[j].Dof; |
|
| if(Dof_P->Type != DOF_UNKNOWN){ | | if(Dof_P->Type != DOF_UNKNOWN) { |
| skip = false; | | skip = false; |
| break; | | break; |
| } | | } |
| } | | } |
| if(skip) return; | | if(skip) return; |
| } | | } |
| } | | } |
| | |
| /* ---------------------------------------------------------------------- */ | | /* ---------------------------------------------------------------------- */ |
| /* G e t F u n c t i o n V a l u e f o r t e s t f u n c t i o n s */ | | /* G e t F u n c t i o n V a l u e f o r t e s t f u n c t i o n s */ |
| /* ---------------------------------------------------------------------- */ | | /* ---------------------------------------------------------------------- */ |
| | |
|
| if (!(Nbr_Equ = QuantityStorageEqu_P->NbrElementaryBasisFunction)) { | | if(!(Nbr_Equ = QuantityStorageEqu_P->NbrElementaryBasisFunction)) { return; } |
| return ; | | |
| } | | |
| | |
| if(Nbr_Equ > NBR_MAX_BASISFUNCTIONS) | | if(Nbr_Equ > NBR_MAX_BASISFUNCTIONS) |
|
| Message::Fatal("Number of basis functions (%d) exceeds NBR_MAX_BASISFUNCTION | | Message::Fatal( |
| S: " | | "Number of basis functions (%d) exceeds NBR_MAX_BASISFUNCTIONS: " |
| "please recompile after changing Interface/ProData.h", Nbr_Eq | | "please recompile after changing src/interface/ProData.h", |
| u); | | Nbr_Equ); |
| | |
| Get_FunctionValue(Nbr_Equ, (void (**)())xFunctionBFEqu, | | Get_FunctionValue(Nbr_Equ, (void (**)())xFunctionBFEqu, |
|
| EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu, | | EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu, |
| QuantityStorageEqu_P, &FI->Type_FormEqu) ; | | QuantityStorageEqu_P, &FI->Type_FormEqu); |
| | |
| /* ---------------------------------------------------------------------- */ | | /* ---------------------------------------------------------------------- */ |
| /* G e t F u n c t i o n V a l u e f o r s h a p e f u n c t i o n s */ | | /* G e t F u n c t i o n V a l u e f o r s h a p e f u n c t i o n s */ |
| /* ---------------------------------------------------------------------- */ | | /* ---------------------------------------------------------------------- */ |
| | |
|
| if (FI->SymmetricalMatrix){ | | if(FI->SymmetricalMatrix) { Nbr_Dof = Nbr_Equ; } |
| Nbr_Dof = Nbr_Equ ; | | else { |
| } | | switch(FI->Type_DefineQuantityDof) { |
| else{ | | case LOCALQUANTITY: |
| switch (FI->Type_DefineQuantityDof) { | | Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction; |
| case LOCALQUANTITY : | | |
| Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ; | | |
| Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof, | | Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof, |
|
| EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof, | | EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof, |
| QuantityStorageDof_P, &FI->Type_FormDof) ; | | QuantityStorageDof_P, &FI->Type_FormDof); |
| break ; | | break; |
| case INTEGRALQUANTITY : | | case INTEGRALQUANTITY: |
| Get_InitElementSource(Element, | | Get_InitElementSource( |
| QuantityStorageDof_P->DefineQuantity->IntegralQuantit | | Element, |
| y.InIndex) ; | | QuantityStorageDof_P->DefineQuantity->IntegralQuantity.InIndex); |
| break ; | | break; |
| case NODOF : | | case NODOF: Nbr_Dof = 1; break; |
| Nbr_Dof = 1 ; | | |
| break ; | | |
| } | | } |
| } | | } |
| | |
| /* ------------------------------------------------------------------- */ | | /* ------------------------------------------------------------------- */ |
| /* G e t I n t e g r a t i o n M e t h o d */ | | /* G e t I n t e g r a t i o n M e t h o d */ |
| /* ------------------------------------------------------------------- */ | | /* ------------------------------------------------------------------- */ |
| | |
| IntegrationCase_P = | | IntegrationCase_P = |
| Get_IntegrationCase(Element, FI->IntegrationCase_L, FI->CriterionIndex); | | Get_IntegrationCase(Element, FI->IntegrationCase_L, FI->CriterionIndex); |
| | |
| /* ------------------------------------------------------------------- */ | | /* ------------------------------------------------------------------- */ |
| /* G e t J a c o b i a n M e t h o d */ | | /* G e t J a c o b i a n M e t h o d */ |
| /* ------------------------------------------------------------------- */ | | /* ------------------------------------------------------------------- */ |
| | |
|
| i = 0 ; | | i = 0; |
| while ((i < FI->NbrJacobianCase) && | | while((i < FI->NbrJacobianCase) && |
| ((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) && | | ((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) && |
| !List_Search | | !List_Search( |
| (((struct Group *)List_Pointer(Problem_S.Group, j)) | | ((struct Group *)List_Pointer(Problem_S.Group, j))->InitialList, |
| ->InitialList, &Element->Region, fcmp_int) ) i++ ; | | &Element->Region, fcmp_int)) |
| | i++; |
| | |
|
| if (i == FI->NbrJacobianCase){ | | if(i == FI->NbrJacobianCase) { |
| Message::Error("Undefined Jacobian in Region %d", Element->Region); | | Message::Error("Undefined Jacobian in Region %d", Element->Region); |
| return; | | return; |
| } | | } |
| | |
|
| Element->JacobianCase = FI->JacobianCase_P0 + i ; | | Element->JacobianCase = FI->JacobianCase_P0 + i; |
| | |
|
| Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*)) | | Get_Jacobian = |
| Get_JacobianFunction(Element->JacobianCase->TypeJacobian, | | (double (*)(struct Element *, MATRIX3x3 *))Get_JacobianFunction( |
| Element->Type, &Type_Dimension) ; | | Element->JacobianCase->TypeJacobian, Element->Type, &Type_Dimension); |
| | |
|
| if (FI->Flag_ChangeCoord) | | if(FI->Flag_ChangeCoord) Get_NodesCoordinatesOfElement(Element); |
| Get_NodesCoordinatesOfElement(Element) ; | | |
| | |
|
| if (Element->JacobianCase->CoefficientIndex < 0){ | | if(Element->JacobianCase->CoefficientIndex < 0) { FI->CoefJac = 1.; } |
| FI->CoefJac = 1.; | | else { |
| } | | Get_ValueOfExpressionByIndex(Element->JacobianCase->CoefficientIndex, NULL, |
| else{ | | 0., 0., 0., &CoefPhys); |
| Get_ValueOfExpressionByIndex(Element->JacobianCase->CoefficientIndex, | | |
| NULL, 0., 0., 0., &CoefPhys) ; | | |
| FI->CoefJac = CoefPhys.Val[0]; | | FI->CoefJac = CoefPhys.Val[0]; |
| } | | } |
| | |
|
| /* ------------------------------------------------------------------------ | | /* ------------------------------------------------------------------------ |
| */ | | */ |
| /* ------------------------------------------------------------------------ | | /* ------------------------------------------------------------------------ |
| */ | | */ |
| /* C o m p u t a t i o n o f E l e m e n t a r y m a t r i x | | /* C o m p u t a t i o n o f E l e m e n t a r y m a t r i x */ |
| */ | | /* ------------------------------------------------------------------------ |
| /* ------------------------------------------------------------------------ | | */ |
| */ | | /* ------------------------------------------------------------------------ |
| /* ------------------------------------------------------------------------ | | */ |
| */ | | |
| | |
| /* Loop on source elements (> 1 only if integral quantity) */ | | /* Loop on source elements (> 1 only if integral quantity) */ |
| | |
|
| while (1) { | | while(1) { |
| | if(FI->Type_DefineQuantityDof == INTEGRALQUANTITY) { |
| | NextElement = Get_NextElementSource(Element->ElementSource); |
| | |
| | if(NextElement) { |
| | /* Get DOF of source element */ |
| | |
| | Get_DofOfElement(Element->ElementSource, |
| | QuantityStorageDof_P->FunctionSpace, |
| | QuantityStorageDof_P, NULL); |
| | |
| | /* Get function value for shape function */ |
| | |
| | Get_NodesCoordinatesOfElement(Element->ElementSource); |
| | Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction; |
| | Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof, |
| | QuantityStorageDof_P->DefineQuantity->IntegralQuantity |
| | .TypeOperatorDof, |
| | QuantityStorageDof_P, |
| | &FI->IntegralQuantityActive.Type_FormDof); |
| | |
|
| if (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) { | | /* Initialize Integral Quantity (integration & jacobian) */ |
| | |
|
| NextElement = Get_NextElementSource(Element->ElementSource) ; | | Cal_InitIntegralQuantity(Element, &FI->IntegralQuantityActive, |
| | QuantityStorageDof_P); |
| if (NextElement) { | | |
| /* Get DOF of source element */ | | |
| | |
| Get_DofOfElement(Element->ElementSource, | | |
| QuantityStorageDof_P->FunctionSpace, | | |
| QuantityStorageDof_P, NULL) ; | | |
| | |
| /* Get function value for shape function */ | | |
| | |
| Get_NodesCoordinatesOfElement(Element->ElementSource) ; | | |
| Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ; | | |
| Get_FunctionValue | | |
| (Nbr_Dof, (void (**)())xFunctionBFDof, | | |
| QuantityStorageDof_P->DefineQuantity->IntegralQuantity.TypeOperatorDof | | |
| , | | |
| QuantityStorageDof_P, &FI->IntegralQuantityActive.Type_FormDof) ; | | |
| | |
| /* Initialize Integral Quantity (integration & jacobian) */ | | |
| | | |
| Cal_InitIntegralQuantity(Element, &FI->IntegralQuantityActive, | | |
| QuantityStorageDof_P); | | |
| } | | } |
| else { | | else { |
|
| break ; | | break; |
| } /* if NextElement */ | | } /* if NextElement */ |
| } /* if INTEGRALQUANTITY */ | | } /* if INTEGRALQUANTITY */ |
| | |
|
| if (FI->SymmetricalMatrix) | | if(FI->SymmetricalMatrix) |
| for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j <= i ; j++) | | for(i = 0; i < Nbr_Equ; i++) |
| for (k = 0 ; k < Current.NbrHar ; k++) Ek[i][j][k] = 0. ; | | for(j = 0; j <= i; j++) |
| | for(k = 0; k < Current.NbrHar; k++) Ek[i][j][k] = 0.; |
| else | | else |
|
| for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j < Nbr_Dof ; j++) | | for(i = 0; i < Nbr_Equ; i++) |
| for (k = 0 ; k < Current.NbrHar ; k++) Ek[i][j][k] = 0. ; | | for(j = 0; j < Nbr_Dof; j++) |
| | for(k = 0; k < Current.NbrHar; k++) Ek[i][j][k] = 0.; |
| switch (IntegrationCase_P->Type) { | | |
| | |
|
| | switch(IntegrationCase_P->Type) { |
| /* ------------------------------------- */ | | /* ------------------------------------- */ |
| /* Q U A D R A T U R E */ | | /* Q U A D R A T U R E */ |
| /* ------------------------------------- */ | | /* ------------------------------------- */ |
| | |
|
| case GAUSS : | | case GAUSS: |
| case GAUSSLEGENDRE : | | case GAUSSLEGENDRE: |
| | |
|
| Quadrature_P = (struct Quadrature*) | | Quadrature_P = (struct Quadrature *)List_PQuery(IntegrationCase_P->Case, |
| List_PQuery(IntegrationCase_P->Case, &Element->Type, fcmp_int); | | &Element->Type, fcmp_int); |
| | |
| if(!Quadrature_P) | | if(!Quadrature_P) |
|
| Message::Error | | Message::Error( |
| ("Unknown type of Element (%s) for Integration method (%s)", | | "Unknown type of Element (%s) for Integration method (%s)", |
| Get_StringForDefine(Element_Type, Element->Type), | | Get_StringForDefine(Element_Type, Element->Type), |
| ((struct IntegrationMethod *) | | ((struct IntegrationMethod *)List_Pointer( |
| List_Pointer(Problem_S.IntegrationMethod, | | Problem_S.IntegrationMethod, |
| EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))- | | EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex)) |
| >Name); | | ->Name); |
| | | |
| Nbr_IntPoints = Quadrature_P->NumberOfPoints ; | | Nbr_IntPoints = Quadrature_P->NumberOfPoints; |
| Get_IntPoint = (void(*)(int,int,double*,double*,double*,double*)) | | Get_IntPoint = (void (*)(int, int, double *, double *, double *, |
| Quadrature_P->Function ; | | double *))Quadrature_P->Function; |
| | |
| for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) { | | |
| | |
|
| | for(i_IntPoint = 0; i_IntPoint < Nbr_IntPoints; i_IntPoint++) { |
| Current.QuadraturePointIndex = i_IntPoint; | | Current.QuadraturePointIndex = i_IntPoint; |
| | |
|
| Get_IntPoint(Nbr_IntPoints, i_IntPoint, | | Get_IntPoint(Nbr_IntPoints, i_IntPoint, &Current.u, &Current.v, |
| &Current.u, &Current.v, &Current.w, &weight) ; | | &Current.w, &weight); |
| | |
|
| if (FI->Flag_ChangeCoord) { | | if(FI->Flag_ChangeCoord) { |
| Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ; | | Get_BFGeoElement(Element, Current.u, Current.v, Current.w); |
| | |
|
| Element->DetJac = Get_Jacobian(Element, &Element->Jac) ; | | Element->DetJac = Get_Jacobian(Element, &Element->Jac); |
| | |
| | if(FI->Flag_InvJac) |
| | Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac, |
| | &Element->Jac, &Element->InvJac); |
| | } |
| | |
|
| if (FI->Flag_InvJac) | | /* Test Functions */ |
| Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac, | | |
| &Element->Jac, &Element->InvJac) ; | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ == |
| } | | CWQ_EXP_TIME_DOF) |
| | | Get_ValueOfExpressionByIndex( |
| /* Test Functions */ | | EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical_Equ, |
| | | QuantityStorage_P0, Current.u, Current.v, Current.w, |
| if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ == CWQ_ | | &CanonicExpression_Equ); |
| EXP_TIME_DOF) | | |
| Get_ValueOfExpressionByIndex | | for(i = 0; i < Nbr_Equ; i++) { |
| (EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical_Equ, | | xFunctionBFEqu[i]( |
| QuantityStorage_P0, Current.u, Current.v, Current.w, &CanonicExpress | | Element, |
| ion_Equ) ; | | QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement + 1, |
| | Current.u, Current.v, Current.w, vBFuEqu[i]); |
| for (i = 0 ; i < Nbr_Equ ; i++) { | | ((void (*)(struct Element *, double *, |
| xFunctionBFEqu[i] | | double *))FI->xChangeOfCoordinatesEqu)(Element, vBFuEqu[i], |
| (Element, | | vBFxEqu[i]); |
| QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement+1, | | |
| Current.u, Current.v, Current.w, vBFuEqu[i]) ; | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != |
| ((void (*)(struct Element*, double*, double*)) | | CWQ_NONE) { |
| FI->xChangeOfCoordinatesEqu) (Element, vBFuEqu[i], vBFxEqu[i]) ; | | V1.Type = Get_ValueFromForm(FI->Type_FormEqu); |
| | | V1.Val[0] = vBFxEqu[i][0]; |
| if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CW | | V1.Val[1] = vBFxEqu[i][1]; |
| Q_NONE){ | | V1.Val[2] = vBFxEqu[i][2]; |
| V1.Type = Get_ValueFromForm(FI->Type_FormEqu); | | V1.Val[MAX_DIM] = 0; |
| V1.Val[0] = vBFxEqu[i][0] ; | | V1.Val[MAX_DIM + 1] = 0; |
| V1.Val[1] = vBFxEqu[i][1] ; | | V1.Val[MAX_DIM + 2] = 0; |
| V1.Val[2] = vBFxEqu[i][2] ; | | |
| V1.Val[MAX_DIM] = 0; | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ == |
| V1.Val[MAX_DIM+1] = 0; | | CWQ_EXP_TIME_DOF) { |
| V1.Val[MAX_DIM+2] = 0; | | switch(EquationTerm_P->Case.LocalTerm.Term |
| | | .OperatorTypeForCanonical_Equ) { |
| if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ == | | case OP_TIME: |
| CWQ_EXP_TIME_DOF){ | | Cal_ProductValue(&CanonicExpression_Equ, &V1, &V2); |
| switch(EquationTerm_P->Case.LocalTerm.Term.OperatorTypeForCanonical | | break; |
| _Equ){ | | case OP_CROSSPRODUCT: |
| case OP_TIME : | | Cal_CrossProductValue(&CanonicExpression_Equ, &V1, &V2); |
| Cal_ProductValue (&CanonicExpression_Equ,&V1,&V2); | | |
| break; | | |
| case OP_CROSSPRODUCT : | | |
| Cal_CrossProductValue (&CanonicExpression_Equ,&V1,&V2); | | |
| break; | | |
| default : | | |
| Message::Error("Unknown operation in Equation"); | | |
| break; | | break; |
|
| } | | default: Message::Error("Unknown operation in Equation"); break; |
| } | | } |
| else if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Eq | | } |
| u == | | else if(EquationTerm_P->Case.LocalTerm.Term |
| CWQ_FCT_DOF){ | | .CanonicalWholeQuantity_Equ == CWQ_FCT_DOF) { |
| ((void(*)(struct Function*, struct Value*, struct Value*)) | | ((void (*)(struct Function *, struct Value *, struct Value *)) |
| EquationTerm_P->Case.LocalTerm.Term.BuiltInFunction_Equ) | | EquationTerm_P->Case.LocalTerm.Term.BuiltInFunction_Equ)( |
| (NULL, &V1, &V2) ; | | NULL, &V1, &V2); |
| } | | } |
| vBFxEqu[i][0] = V2.Val[0]; | | vBFxEqu[i][0] = V2.Val[0]; |
| vBFxEqu[i][1] = V2.Val[1]; | | vBFxEqu[i][1] = V2.Val[1]; |
| vBFxEqu[i][2] = V2.Val[2]; | | vBFxEqu[i][2] = V2.Val[2]; |
| } | | } |
| | | |
| } /* for Nbr_Equ */ | | } /* for Nbr_Equ */ |
| | | |
| /* Shape Functions (+ surrounding expression) */ | | /* Shape Functions (+ surrounding expression) */ |
| | | |
| Current.Element = Element ; | | Current.Element = Element; |
| Cal_vBFxDof(EquationTerm_P, FI, | | if(Cal_vBFxDof(EquationTerm_P, FI, QuantityStorage_P0, |
| QuantityStorage_P0, QuantityStorageDof_P, | | QuantityStorageDof_P, Nbr_Dof, xFunctionBFDof, vBFxEqu, |
| Nbr_Dof, xFunctionBFDof, vBFxEqu, vBFxDof); | | vBFxDof)) { |
| | | Factor = |
| Factor = FI->CoefJac * | | FI->CoefJac * |
| ((FI->Flag_ChangeCoord) ? weight * fabs(Element->DetJac) : weight) ; | | ((FI->Flag_ChangeCoord) ? weight * fabs(Element->DetJac) : weight); |
| | | |
| /* Product and assembly in elementary submatrix (k?-1.:1.)* | | /* Product and assembly in elementary submatrix (k?-1.:1.)* */ |
| */ | | if(FI->SymmetricalMatrix) |
| if (FI->SymmetricalMatrix) | | for(i = 0; i < Nbr_Equ; i++) |
| for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j <= i ; j++) | | for(j = 0; j <= i; j++) |
| for (k = 0 ; k < Current.NbrHar ; k++) | | for(k = 0; k < Current.NbrHar; k++) |
| Ek[i][j][k] += Factor * | | Ek[i][j][k] += |
| ((double (*)(double*, double*)) | | Factor * ((double (*)(double *, double *))FI->Cal_Productx)( |
| FI->Cal_Productx) (vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM*k])) ; | | vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM * k])); |
| else | | else |
| for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j < Nbr_Dof ; j++) | | for(i = 0; i < Nbr_Equ; i++) |
| for (k = 0 ; k < Current.NbrHar ; k++) | | for(j = 0; j < Nbr_Dof; j++) |
| Ek[i][j][k] += Factor * | | for(k = 0; k < Current.NbrHar; k++) |
| ((double (*)(double*, double*)) | | Ek[i][j][k] += |
| FI->Cal_Productx) (vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM*k])); | | Factor * ((double (*)(double *, double *))FI->Cal_Productx)( |
| | vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM * k])); |
| | } |
| | |
|
| } /* for i_IntPoint ... */ | | } /* for i_IntPoint ... */ |
| break ; /* case GAUSS */ | | break; /* case GAUSS */ |
| | |
| /* ------------------------------------- */ | | /* ------------------------------------- */ |
| /* A N A L Y T I C */ | | /* A N A L Y T I C */ |
| /* ------------------------------------- */ | | /* ------------------------------------- */ |
| | |
|
| case ANALYTIC : | | case ANALYTIC: |
| | |
|
| if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == |
| CWQ_DOF) { | | CWQ_DOF) { |
| Factor = 1. ; | | Factor = 1.; |
| } | | } |
| if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == | | if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity == |
| CWQ_EXP_TIME_DOF) { | | CWQ_EXP_TIME_DOF) { |
| if (EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical >= 0) | | if(EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical >= |
| { | | 0) { |
| Get_ValueOfExpression | | Get_ValueOfExpression( |
| ((struct Expression *)List_Pointer | | (struct Expression *)List_Pointer( |
| (Problem_S.Expression, | | Problem_S.Expression, |
| EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical), | | EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical), |
| QuantityStorage_P0, 0., 0., 0., &CoefPhys) ; | | QuantityStorage_P0, 0., 0., 0., &CoefPhys); |
| Factor = CoefPhys.Val[0] ; | | Factor = CoefPhys.Val[0]; |
| } | | } |
| } | | } |
| else { | | else { |
|
| Message::Error("Bad expression for full analytic integration"); | | Message::Error("Bad expression for full analytic integration"); |
| } | | } |
| | |
|
| if (FI->SymmetricalMatrix) { | | if(FI->SymmetricalMatrix) { |
| for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j <= i ; j++) | | for(i = 0; i < Nbr_Equ; i++) |
| Ek[i][j][0] = Factor * | | for(j = 0; j <= i; j++) |
| Cal_AnalyticIntegration | | Ek[i][j][0] = |
| (Element, (void (*)())xFunctionBFEqu[i], (void (*)())xFunctionBFEqu[j | | Factor * Cal_AnalyticIntegration( |
| ], i, j, | | Element, (void (*)())xFunctionBFEqu[i], |
| FI->Cal_Productx) ; | | (void (*)())xFunctionBFEqu[j], i, j, FI->Cal_Productx); |
| } | | } |
| else { | | else { |
|
| switch (FI->Type_DefineQuantityDof) { | | switch(FI->Type_DefineQuantityDof) { |
| case LOCALQUANTITY : | | case LOCALQUANTITY: |
| for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j < Nbr_Dof ; j++) | | for(i = 0; i < Nbr_Equ; i++) |
| Ek[i][j][0] = Factor * | | for(j = 0; j < Nbr_Dof; j++) |
| Cal_AnalyticIntegration | | Ek[i][j][0] = Factor * Cal_AnalyticIntegration( |
| (Element, (void (*)())xFunctionBFEqu[i], (void (*)())xFunctionBFDof | | Element, (void (*)())xFunctionBFEqu[i], |
| [j], i, j, | | (void (*)())xFunctionBFDof[j], i, j, |
| FI->Cal_Productx) ; | | FI->Cal_Productx); |
| break; | | break; |
| default : | | default: |
| Message::Error("Exterior analytical integration not implemented"); | | Message::Error("Exterior analytical integration not implemented"); |
| break; | | break; |
| } | | } |
| } | | } |
| | | |
| break ; /* case ANALYTIC */ | | break; /* case ANALYTIC */ |
| | | |
| default : | | |
| Message::Error | | |
| ("Unknown type of Integration method (%s)", | | |
| ((struct IntegrationMethod *) | | |
| List_Pointer(Problem_S.IntegrationMethod, | | |
| EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->N | | |
| ame); | | |
| break; | | |
| | |
|
| | default: |
| | Message::Error("Unknown type of Integration method (%s)", |
| | ((struct IntegrationMethod *)List_Pointer( |
| | Problem_S.IntegrationMethod, |
| | EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex)) |
| | ->Name); |
| | break; |
| } | | } |
| | |
| /* Complete elementary matrix if symmetrical */ | | /* Complete elementary matrix if symmetrical */ |
| /* ----------------------------------------- */ | | /* ----------------------------------------- */ |
| | |
|
| if (FI->SymmetricalMatrix) | | if(FI->SymmetricalMatrix) |
| for (i = 1 ; i < Nbr_Equ ; i++) | | for(i = 1; i < Nbr_Equ; i++) |
| for (j = 0 ; j < i ; j++) | | for(j = 0; j < i; j++) |
| for (k = 0 ; k < Current.NbrHar ; k++) | | for(k = 0; k < Current.NbrHar; k++) Ek[j][i][k] = Ek[i][j][k]; |
| Ek[j][i][k] = Ek[i][j][k] ; | | |
| | |
| if(Message::GetVerbosity() == 10) { | | if(Message::GetVerbosity() == 10) { |
|
| printf("Galerkin = ") ; | | printf("Galerkin = "); |
| for (j = 0 ; j < Nbr_Dof ; j++) | | for(j = 0; j < Nbr_Dof; j++) |
| Print_DofNumber(QuantityStorageDof_P->BasisFunction[j].Dof) ; | | Print_DofNumber(QuantityStorageDof_P->BasisFunction[j].Dof); |
| printf("\n") ; | | printf("\n"); |
| for (i = 0 ; i < Nbr_Equ ; i++) { | | for(i = 0; i < Nbr_Equ; i++) { |
| Print_DofNumber(QuantityStorageEqu_P->BasisFunction[i].Dof) ; | | Print_DofNumber(QuantityStorageEqu_P->BasisFunction[i].Dof); |
| printf("[ ") ; | | printf("[ "); |
| for (j = 0 ; j < Nbr_Dof ; j++) { | | for(j = 0; j < Nbr_Dof; j++) { |
| printf("(") ; | | printf("("); |
| for(k = 0 ; k < Current.NbrHar ; k++) printf("% .5e ", Ek[i][j][k]) ; | | for(k = 0; k < Current.NbrHar; k++) printf("% .5e ", Ek[i][j][k]); |
| printf(")") ; | | printf(")"); |
| } | | } |
| printf("]\n") ; | | printf("]\n"); |
| } | | } |
| } | | } |
| | |
| // store unassembled RHS in DofData, per element | | // store unassembled RHS in DofData, per element |
| #if 0 | | #if 0 |
| for (j = 0 ; j < Nbr_Dof ; j++){ | | for (j = 0 ; j < Nbr_Dof ; j++){ |
| if(QuantityStorageDof_P->BasisFunction[j].Dof->Type == DOF_FIXED && | | if(QuantityStorageDof_P->BasisFunction[j].Dof->Type == DOF_FIXED && |
| QuantityStorageDof_P->BasisFunction[j].Dof->Entity == 0){ | | QuantityStorageDof_P->BasisFunction[j].Dof->Entity == 0){ |
| | |
| std::vector<std::pair<int, double> > &vec = | | std::vector<std::pair<int, double> > &vec = |
| | |
| skipping to change at line 936 | | skipping to change at line 969 |
|---|
| } | | } |
| } | | } |
| } | | } |
| printf("\n") ; | | printf("\n") ; |
| } | | } |
| } | | } |
| #endif | | #endif |
| | |
| /* Assembly in global matrix */ | | /* Assembly in global matrix */ |
| /* ------------------------- */ | | /* ------------------------- */ |
|
| if (!Current.flagAssDiag) /*+++prov*/ | | if(!Current.flagAssDiag) /*+++prov*/ |
| for (i = 0 ; i < Nbr_Equ ; i++) | | for(i = 0; i < Nbr_Equ; i++) |
| for (j = 0 ; j < Nbr_Dof ; j++) | | for(j = 0; j < Nbr_Dof; j++) |
| ((void (*)(struct Dof*, struct Dof*, double*)) | | ((void (*)(struct Dof *, struct Dof *, |
| FI->Function_AssembleTerm) | | double *))FI->Function_AssembleTerm)( |
| (QuantityStorageEqu_P->BasisFunction[i].Dof, | | QuantityStorageEqu_P->BasisFunction[i].Dof, |
| QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]) ; | | QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]); |
| else if (Current.flagAssDiag == 1) { | | else if(Current.flagAssDiag == 1) { |
| for (i = 0 ; i < Nbr_Equ ; i++) { | | for(i = 0; i < Nbr_Equ; i++) { |
| /* for (j = 0 ; j < Nbr_Dof ; j++)*/ | | /* for (j = 0 ; j < Nbr_Dof ; j++)*/ |
| j = i; | | j = i; |
| ((void (*)(struct Dof*, struct Dof*, double*)) | | ((void (*)(struct Dof *, struct Dof *, |
| FI->Function_AssembleTerm) | | double *))FI->Function_AssembleTerm)( |
| (QuantityStorageEqu_P->BasisFunction[i].Dof, | | QuantityStorageEqu_P->BasisFunction[i].Dof, |
| QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]) ; | | QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]); |
| } | | } |
| } | | } |
| else if (Current.flagAssDiag == 2) { | | else if(Current.flagAssDiag == 2) { |
| for (i = 0 ; i < Nbr_Equ ; i++) { | | for(i = 0; i < Nbr_Equ; i++) { |
| /* for (j = 0 ; j < Nbr_Dof ; j++)*/ | | /* for (j = 0 ; j < Nbr_Dof ; j++)*/ |
| j = i; | | j = i; |
| if (QuantityStorageEqu_P->BasisFunction[i].Dof->Type == DOF_UNKNOWN | | if(QuantityStorageEqu_P->BasisFunction[i].Dof->Type == DOF_UNKNOWN && |
| && | | assDiag_done.find( |
| assDiag_done.find | | QuantityStorageEqu_P->BasisFunction[i].Dof->Case.Unknown.NumDof - |
| (QuantityStorageEqu_P->BasisFunction[i].Dof->Case.Unknown.NumDof-1) | | 1) == assDiag_done.end()) { |
| == assDiag_done.end()) { | | assDiag_done[QuantityStorageEqu_P->BasisFunction[i] |
| assDiag_done | | .Dof->Case.Unknown.NumDof - |
| [QuantityStorageEqu_P->BasisFunction[i].Dof->Case.Unknown.NumDof-1] | | 1] = true; |
| = true; | | |
| Ek[i][j][0] = 1.; | | Ek[i][j][0] = 1.; |
|
| for (k = 1 ; k < Current.NbrHar ; k++) Ek[i][j][k] = 0. ; | | for(k = 1; k < Current.NbrHar; k++) Ek[i][j][k] = 0.; |
| ((void (*)(struct Dof*, struct Dof*, double*)) | | ((void (*)(struct Dof *, struct Dof *, |
| FI->Function_AssembleTerm) | | double *))FI->Function_AssembleTerm)( |
| (QuantityStorageEqu_P->BasisFunction[i].Dof, | | QuantityStorageEqu_P->BasisFunction[i].Dof, |
| QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]) ; | | QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]); |
| } | | } |
| } | | } |
| } | | } |
| | |
| /* Exit while(1) directly if not integral quantity */ | | /* Exit while(1) directly if not integral quantity */ |
| | |
|
| if (FI->Type_DefineQuantityDof != INTEGRALQUANTITY) break ; | | if(FI->Type_DefineQuantityDof != INTEGRALQUANTITY) break; |
| | |
|
| } /* while (1) ... */ | | } /* while (1) ... */ |
| | |
| Current.flagAssDiag = 0; /*+++prov*/ | | Current.flagAssDiag = 0; /*+++prov*/ |
|
| | | |
| } | | } |
| | | |
| End of changes. 118 change blocks. |
|---|
| 670 lines changed or deleted | | 639 lines changed or added |
|---|
"Fossies" - the Fresh Open Source Software Archive 