#define RCSID "$Id: Cal_GalerkinTermOfFemEquation.c,v 1.25 2006/02/26 00:42:54 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>.
*
* Contributor(s):
* Johan Gyselinck
* Ruth Sabariego
*/
#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "Get_DofOfElement.h"
#include "Get_Geometry.h"
#include "Data_DefineE.h"
#include "GeoData.h"
#include "CurrentData.h"
#include "Tools.h"
#include "F_FMM.h"
void Cal_InitGalerkinTermOfFemEquation_MHJacNL(struct EquationTerm * EquationTerm_P) ;
void Cal_GalerkinTermOfFemEquation_MHJacNL(struct Element * Element,
struct EquationTerm * EquationTerm_P,
struct QuantityStorage * QuantityStorage_P0) ;
/* ------------------------------------------------------------------------ */
/* 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,
struct QuantityStorage * QuantityStorage_P0,
struct QuantityStorage * QuantityStorageNoDof,
struct Dof * DofForNoDof_P) {
struct FemLocalTermActive * FI ;
extern int MH_Moving_Matrix_simple, MH_Moving_Matrix_probe, MH_Moving_Matrix_separate;
GetDP_Begin("Cal_InitGalerkinTermOfFemEquation");
FI = EquationTerm_P->Case.LocalTerm.Active ;
FI->QuantityStorageEqu_P = QuantityStorage_P0 +
EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ;
Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu,
FI->QuantityStorageEqu_P, &FI->Type_FormEqu) ;
if (EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof >= 0) {
FI->QuantityStorageDof_P = QuantityStorage_P0 +
EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof ;
FI->Type_DefineQuantityDof = FI->QuantityStorageDof_P->DefineQuantity->Type ;
}
else {
FI->QuantityStorageDof_P = QuantityStorageNoDof ;
FI->Type_DefineQuantityDof = NODOF ;
FI->DofForNoDof_P = DofForNoDof_P ;
Dof_InitDofForNoDof(DofForNoDof_P, Current.NbrHar) ;
QuantityStorageNoDof->BasisFunction[0].Dof = DofForNoDof_P ;
}
/* Warning: not correct if nonsymmetrical tensor in expression */
FI->SymmetricalMatrix =
(EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ==
EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexDof) &&
(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu ==
EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof) ;
if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE)
FI->SymmetricalMatrix = 0 ;
if (FI->SymmetricalMatrix) {
FI->Type_FormDof = FI->Type_FormEqu ;
}
else {
switch (FI->Type_DefineQuantityDof) {
case LOCALQUANTITY :
Get_InitFunctionValue(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
FI->QuantityStorageDof_P, &FI->Type_FormDof) ;
break ;
case INTEGRALQUANTITY :
if(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof != NOOP){
Msg(GERROR, "No operator can act on an Integral Quantity");
}
FI->Type_FormDof = VECTOR ; /* we don't know the type a priori */
FI->IntegralQuantityActive.IntegrationCase_L =
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
FI->QuantityStorageDof_P->DefineQuantity->
IntegralQuantity.IntegrationMethodIndex)) ->IntegrationCase ;
FI->IntegralQuantityActive.CriterionIndex =
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
FI->QuantityStorageDof_P->DefineQuantity->
IntegralQuantity.IntegrationMethodIndex)) ->CriterionIndex ;
FI->IntegralQuantityActive.JacobianCase_L =
((struct JacobianMethod *)
List_Pointer(Problem_S.JacobianMethod,
FI->QuantityStorageDof_P->DefineQuantity->
IntegralQuantity.JacobianMethodIndex)) ->JacobianCase ;
break ;
case NODOF :
FI->Type_FormDof = SCALAR ;
break ;
}
}
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 */
/* -------------------------------------------- */
if(EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex < 0)
Msg(GERROR, "Integration method missing in equation term");
FI->IntegrationCase_L =
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))
->IntegrationCase ;
FI->CriterionIndex =
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))
->CriterionIndex ;
/* G e t J a c o b i a n M e t h o d */
/* -------------------------------------- */
if(EquationTerm_P->Case.LocalTerm.JacobianMethodIndex < 0)
Msg(GERROR, "Jacobian method missing in equation term");
FI->JacobianCase_L =
((struct JacobianMethod *)
List_Pointer(Problem_S.JacobianMethod,
EquationTerm_P->Case.LocalTerm.JacobianMethodIndex))
->JacobianCase ;
FI->JacobianCase_P0 =
(FI->NbrJacobianCase = List_Nbr(FI->JacobianCase_L)) ?
(struct JacobianCase*)List_Pointer(FI->JacobianCase_L, 0) : NULL ;
FI->Flag_ChangeCoord =
( FI->SymmetricalMatrix ||
!(
( (FI->Type_FormEqu == FORM0 || FI->Type_FormEqu == FORM0P) &&
(FI->Type_FormDof == FORM3 || FI->Type_FormDof == FORM3P) ) ||
/*
( (FI->Type_FormEqu == FORM1 || FI->Type_FormEqu == FORM1P) &&
(FI->Type_FormDof == FORM2 || FI->Type_FormDof == FORM2P) ) ||
( (FI->Type_FormEqu == FORM2 || FI->Type_FormEqu == FORM2P) &&
(FI->Type_FormDof == FORM1 || FI->Type_FormDof == FORM1P) ) ||
*/
( (FI->Type_FormEqu == FORM3 || FI->Type_FormEqu == FORM3P) &&
(FI->Type_FormDof == FORM0 || FI->Type_FormDof == FORM0P) )
)
)
|| /* For operators on VECTOR's (To be extended) */
(FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR)
||
(FI->Type_DefineQuantityDof == INTEGRALQUANTITY) ;
if (FI->Flag_ChangeCoord){
FI->Flag_InvJac = ( (FI->Type_FormEqu == FORM1) ||
(!FI->SymmetricalMatrix && (FI->Type_FormDof == FORM1)) ||
/* For operators on VECTOR's (To be extended) */
(FI->Type_FormEqu == VECTOR || FI->Type_FormDof == VECTOR) ||
(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 */
/* ---------------------------------------------- */
FI->xChangeOfCoordinatesEqu =
(void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormEqu) ;
if (!FI->SymmetricalMatrix)
FI->xChangeOfCoordinatesDof =
(void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormDof) ;
else
FI->xChangeOfCoordinatesDof =
(void (*)())Get_ChangeOfCoordinates(0, FI->Type_FormDof) ; /* Used only for transfer */
/* G e t C a l _ P r o d u c t x */
/* -------------------------------- */
switch (FI->Type_FormEqu) {
case FORM1 : case FORM1S :
case FORM2 : case FORM2P : case FORM2S :
case VECTOR :
FI->Cal_Productx = (double (*)())Cal_Product123 ; break ;
case FORM1P :
case VECTORP :
FI->Cal_Productx = (double (*)())Cal_Product3 ; break ;
case FORM0 :
case FORM3 : case FORM3P :
case SCALAR :
FI->Cal_Productx = (double (*)())Cal_Product1 ; break ;
default :
Msg(GERROR, "Unknown type of Form (%d)", FI->Type_FormEqu);
}
/* 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) {
case NODT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NoDt ; break;
case DT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_Dt ; break;
case DTDOF_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDof ; break;
case DTDT_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDt ; break;
case DTDTDOF_: FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtDtDof; break;
case JACNL_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_JacNL ; break;
case NEVERDT_: FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_NeverDt; break;
case DTNL_ : FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_DtNL ; break;
default : Msg(GERROR, "Unknown type of Operator for Galerkin term (%d)",
EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative);
}
if (MH_Moving_Matrix_simple) {
/* Msg (INFO, "AssembleTerm_MH_Moving") ; */
FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_simple ;
}
if (MH_Moving_Matrix_probe) {
/* Msg (INFO, "AssembleTerm_MH_Moving") ; */
FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_probe ;
}
if (MH_Moving_Matrix_separate) {
/* Msg (INFO, "AssembleTerm_MH_Moving") ; */
FI->Function_AssembleTerm = (void (*)())Cal_AssembleTerm_MH_Moving_separate ;
}
/* initialisation of MHJacNL-term (nonlinear multi-harmonics) if necessary */
Cal_InitGalerkinTermOfFemEquation_MHJacNL(EquationTerm_P);
/* Full_Matrix */
if (EquationTerm_P->Case.LocalTerm.Full_Matrix) {
FI->Full_Matrix = 1;
FI->FirstElements = List_Create(20, 10, sizeof(struct FirstElement)) ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* 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,
struct EquationTerm * EquationTerm_P,
struct QuantityStorage * QuantityStorage_P0) {
struct FemLocalTermActive * FI ;
struct QuantityStorage * QuantityStorageEqu_P, * QuantityStorageDof_P ;
struct IntegrationCase * IntegrationCase_P ;
struct Quadrature * Quadrature_P ;
struct Value vBFxDof [NBR_MAX_BASISFUNCTIONS], CoefPhys ;
struct Value CanonicExpression_Equ, V1, V2;
int Nbr_Equ, Nbr_Dof = 0;
int i, j, k, Type_Dimension, Nbr_IntPoints, i_IntPoint ;
int NextElement ;
double weight, Factor = 1. ;
double vBFuEqu [NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ;
double vBFxEqu [NBR_MAX_BASISFUNCTIONS] [MAX_DIM] ;
double Ek [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_HARMONIC] ;
void (*xFunctionBFEqu[NBR_MAX_BASISFUNCTIONS])
(struct Element * Element, int NumEntity,
double u, double v, double w, double Value[] ) ;
void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
(struct Element * Element, int NumEntity,
double u, double v, double w, double Value[] ) ;
double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
void (*Get_IntPoint)(int,int,double*,double*,double*,double*);
extern int Flag_RHS;
GetDP_Begin("Cal_GalerkinTermOfFemEquation");
FI = EquationTerm_P->Case.LocalTerm.Active ;
/* treatment of MHJacNL-term in separate routine */
if (FI->MHJacNL) {
/* if only the RHS of the system is to be calculated
(in case of adaptive relaxation of the Newton-Raphson scheme)
the (expensive and redundant) calculation of this term can be skipped */
if (!Flag_RHS)
Cal_GalerkinTermOfFemEquation_MHJacNL(Element, EquationTerm_P, QuantityStorage_P0) ;
GetDP_End ;
}
QuantityStorageEqu_P = FI->QuantityStorageEqu_P ;
QuantityStorageDof_P = FI->QuantityStorageDof_P ;
/* ---------------------------------------------------------------------- */
/* 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)) {
GetDP_End ;
}
Get_FunctionValue(Nbr_Equ, (void (**)())xFunctionBFEqu,
EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu,
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 */
/* ---------------------------------------------------------------------- */
if (FI->SymmetricalMatrix){
Nbr_Dof = Nbr_Equ ;
}
else{
switch (FI->Type_DefineQuantityDof) {
case LOCALQUANTITY :
Nbr_Dof = QuantityStorageDof_P->NbrElementaryBasisFunction ;
Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
QuantityStorageDof_P, &FI->Type_FormDof) ;
break ;
case INTEGRALQUANTITY :
Get_InitElementSource(Element,
QuantityStorageDof_P->DefineQuantity->IntegralQuantity.InIndex) ;
break ;
case NODOF :
Nbr_Dof = 1 ;
break ;
}
}
/* ------------------------------------------------------------------- */
/* G e t I n t e g r a t i o n M e t h o d */
/* ------------------------------------------------------------------- */
IntegrationCase_P =
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 */
/* ------------------------------------------------------------------- */
i = 0 ;
while ((i < FI->NbrJacobianCase) &&
((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) &&
!List_Search
(((struct Group *)List_Pointer(Problem_S.Group, j))
->InitialList, &Element->Region, fcmp_int) ) i++ ;
if (i == FI->NbrJacobianCase)
Msg(GERROR, "Undefined Jacobian in Region %d", Element->Region);
Element->JacobianCase = FI->JacobianCase_P0 + i ;
Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
Element->Type, &Type_Dimension) ;
if (FI->Flag_ChangeCoord)
Get_NodesCoordinatesOfElement(Element) ;
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
/* 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) */
while (1) {
if (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) {
if (!Flag_FMM)
NextElement = Get_NextElementSource(Element->ElementSource) ;
else{
NextElement = Get_NextElementSourceNeighbour(Element) ;
}
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 {
break ;
} /* if NextElement */
} /* if INTEGRALQUANTITY */
if (FI->SymmetricalMatrix)
for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j <= i ; j++)
for (k = 0 ; k < Current.NbrHar ; k++) Ek[i][j][k] = 0. ;
else
for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j < Nbr_Dof ; j++)
for (k = 0 ; k < Current.NbrHar ; k++) Ek[i][j][k] = 0. ;
switch (IntegrationCase_P->Type) {
/* ------------------------------------- */
/* Q U A D R A T U R E */
/* ------------------------------------- */
case GAUSS :
case GAUSSLEGENDRE :
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)",
Get_StringForDefine(Element_Type, Element->Type),
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
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,
&Current.u, &Current.v, &Current.w, &weight) ;
if (FI->Flag_ChangeCoord) {
Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ;
Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
if (FI->Flag_InvJac)
Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
&Element->Jac, &Element->InvJac) ;
}
/* Test Functions */
if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE)
Get_ValueOfExpressionByIndex
(EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical_Equ,
QuantityStorage_P0, Current.u, Current.v, Current.w, &CanonicExpression_Equ) ;
for (i = 0 ; i < Nbr_Equ ; i++) {
xFunctionBFEqu[i]
(Element,
QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement+1,
Current.u, Current.v, Current.w, vBFuEqu[i]) ;
((void (*)(struct Element*, double*, double*))
FI->xChangeOfCoordinatesEqu) (Element, vBFuEqu[i], vBFxEqu[i]) ;
if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE){
V1.Type = Get_ValueFromForm(FI->Type_FormEqu);
V1.Val[0] = vBFxEqu[i][0] ;
V1.Val[1] = vBFxEqu[i][1] ;
V1.Val[2] = vBFxEqu[i][2] ;
V1.Val[MAX_DIM] = 0;
V1.Val[MAX_DIM+1] = 0;
V1.Val[MAX_DIM+2] = 0;
switch(EquationTerm_P->Case.LocalTerm.Term.OperatorTypeForCanonical_Equ){
case OP_TIME : Cal_ProductValue (&CanonicExpression_Equ,&V1,&V2); break;
case OP_CROSSPRODUCT : Cal_CrossProductValue (&CanonicExpression_Equ,&V1,&V2); break;
default : Msg(GERROR, "Unknown operation in Equation");
}
vBFxEqu[i][0] = V2.Val[0];
vBFxEqu[i][1] = V2.Val[1];
vBFxEqu[i][2] = V2.Val[2];
}
} /* for Nbr_Equ */
/* Shape Functions (+ surrounding expression) */
Current.Element = Element ;
Cal_vBFxDof(EquationTerm_P, FI,
QuantityStorage_P0, QuantityStorageDof_P,
Nbr_Dof, xFunctionBFDof, vBFxEqu, vBFxDof);
Factor = (FI->Flag_ChangeCoord) ? weight * fabs(Element->DetJac) : weight ;
/* Product and assembly in elementary submatrix (k?-1.:1.)* */
if (FI->SymmetricalMatrix)
for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j <= i ; j++)
for (k = 0 ; k < Current.NbrHar ; k++)
Ek[i][j][k] += Factor *
((double (*)(double*, double*))
FI->Cal_Productx) (vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM*k])) ;
else
for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j < Nbr_Dof ; j++)
for (k = 0 ; k < Current.NbrHar ; k++)
Ek[i][j][k] += Factor *
((double (*)(double*, double*))
FI->Cal_Productx) (vBFxEqu[i], &(vBFxDof[j].Val[MAX_DIM*k])) ;
} /* for i_IntPoint ... */
break ; /* case GAUSS */
/* ------------------------------------- */
/* A N A L Y T I C */
/* ------------------------------------- */
case ANALYTIC :
if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity ==
CWQ_DOF) {
Factor = 1. ;
}
if (EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity ==
CWQ_EXP_TIME_DOF) {
if (EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical >= 0) {
Get_ValueOfExpression
(Problem_Expression0 +
EquationTerm_P->Case.LocalTerm.Term.ExpressionIndexForCanonical,
QuantityStorage_P0, 0., 0., 0., &CoefPhys) ;
Factor = CoefPhys.Val[0] ;
}
}
else {
Msg(GERROR, "Bad expression for full analytic integration");
}
if (FI->SymmetricalMatrix) {
for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j <= i ; j++)
Ek[i][j][0] = Factor *
Cal_AnalyticIntegration
(Element, (void (*)())xFunctionBFEqu[i], (void (*)())xFunctionBFEqu[j], i, j,
FI->Cal_Productx) ;
}
else {
switch (FI->Type_DefineQuantityDof) {
case LOCALQUANTITY :
for (i = 0 ; i < Nbr_Equ ; i++) for (j = 0 ; j < Nbr_Dof ; j++)
Ek[i][j][0] = Factor *
Cal_AnalyticIntegration
(Element, (void (*)())xFunctionBFEqu[i], (void (*)())xFunctionBFDof[j], i, j,
FI->Cal_Productx) ;
break;
default :
Msg(GERROR, "Exterior analytical integration not implemented");
break;
}
}
break ; /* case ANALYTIC */
default :
Msg(GERROR, "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 */
/* ----------------------------------------- */
if (FI->SymmetricalMatrix)
for (i = 1 ; i < Nbr_Equ ; i++)
for (j = 0 ; j < i ; j++)
for (k = 0 ; k < Current.NbrHar ; k++)
Ek[j][i][k] = Ek[i][j][k] ;
if(Flag_VERBOSE == 10) {
printf("Galerkin = ") ;
for (j = 0 ; j < Nbr_Dof ; j++)
Print_DofNumber(QuantityStorageDof_P->BasisFunction[j].Dof) ;
printf("\n") ;
for (i = 0 ; i < Nbr_Equ ; i++) {
Print_DofNumber(QuantityStorageEqu_P->BasisFunction[i].Dof) ;
printf("[ ") ;
for (j = 0 ; j < Nbr_Dof ; j++) {
printf("(") ;
for(k = 0 ; k < Current.NbrHar ; k++) printf("% .5e ", Ek[i][j][k]) ;
printf(")") ;
}
printf("]\n") ;
}
}
/* Assembly in global matrix */
/* ------------------------- */
for (i = 0 ; i < Nbr_Equ ; i++)
for (j = 0 ; j < Nbr_Dof ; j++)
((void (*)(struct Dof*, struct Dof*, double*))
FI->Function_AssembleTerm)
(QuantityStorageEqu_P->BasisFunction[i].Dof,
QuantityStorageDof_P->BasisFunction[j].Dof, Ek[i][j]) ;
/* Exit while(1) directly if not integral quantity */
if (FI->Type_DefineQuantityDof != INTEGRALQUANTITY) break ;
} /* while (1) ... */
GetDP_End ;
}
syntax highlighted by Code2HTML, v. 0.9.1