#define RCSID "$Id: Cal_deRhamTermOfFemEquation.c,v 1.20 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>.
*/
#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "Get_DofOfElement.h"
#include "Get_Geometry.h"
#include "Get_Cells.h"
#include "Pos_Search.h" /* xyz2uvw */
#include "Numeric.h"
#include "Data_DefineE.h"
#include "GeoData.h"
#include "CurrentData.h"
#include "Tools.h"
/* ------------------------------------------------------------------------ */
/* C a l _ I n i t D e R h a m T e r m O f F e m E q u a t i o n */
/* ------------------------------------------------------------------------ */
void Cal_InitdeRhamTermOfFemEquation(struct EquationTerm * EquationTerm_P,
struct QuantityStorage * QuantityStorage_P0,
struct QuantityStorage * QuantityStorageNoDof,
struct Dof * DofForNoDof_P) {
struct FemLocalTermActive * FI ;
GetDP_Begin("Cal_InitdeRhamTermOfFemEquation");
FI = EquationTerm_P->Case.LocalTerm.Active ;
FI->QuantityStorageEqu_P = QuantityStorage_P0 +
EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ;
if(EquationTerm_P->Case.LocalTerm.Term.TypeOperatorEqu != NOOP)
Msg(GERROR, "An aperator cannot act on the \"test function\" in a de Rham Map");
FI->Type_FormEqu = -1 ;
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 ;
}
FI->SymmetricalMatrix = 0 ;
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, "A differential operator cannot 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 */
/* -------------------------------------------- */
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 */
/* -------------------------------------- */
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 = 1 ;
FI->Flag_InvJac =
(FI->Type_FormDof == FORM1) ||
(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 = NULL ;
FI->xChangeOfCoordinatesDof =
(void (*)())Get_ChangeOfCoordinates(FI->Flag_ChangeCoord, FI->Type_FormDof) ;
FI->Cal_Productx = NULL ;
/* 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;
default : Msg(GERROR, "Unknown type of operator for de Rham (%d)",
EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative);
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* C a l _ D e R h a m T e r m O f F e m E q u a t i o n */
/* ------------------------------------------------------------------------ */
void Cal_deRhamTermOfFemEquation(struct Element * Element,
struct EquationTerm * EquationTerm_P,
struct QuantityStorage * QuantityStorage_P0) {
struct Element Cells[NBR_MAX_DERHAM_CELLS];
struct FemLocalTermActive * FI ;
struct QuantityStorage * QuantityStorageEqu_P, * QuantityStorageDof_P ;
struct IntegrationCase * IntegrationCase_P ;
struct Quadrature * Quadrature_P ;
struct Value vBFxDof[NBR_MAX_BASISFUNCTIONS] ;
struct Value Cell_Vector[NBR_MAX_DERHAM_CELLS];
struct Value Cell_Value[NBR_MAX_BASISFUNCTIONS][NBR_MAX_DERHAM_CELLS];
struct Group * Group_P ;
double u, v, w, weight, Factor ;
double u2, v2, w2, x2, y2, z2 ;
double Ek [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_BASISFUNCTIONS] [NBR_MAX_HARMONIC] ;
int Nbr_Equ, Nbr_Dof = 0 ;
int i, j, k, n, ii, Nbr_IntPoints, i_IntPoint, i_Cell ;
int Type_Dimension, Cell_Type_Dimension ;
int Nbr_Cells, Cell_RelativeJacobianType ;
int Nbr_Ent1 = 0, Nbr_Ent2 = 0, *Dk = NULL ;
int Product ;
void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
(struct Element * Element, int NumEntity,
double u, double v, double w, double Value[] ) ;
double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
double (*Cell_Get_Jacobian)(struct Element*, MATRIX3x3*) = 0;
void (*Get_IntPoint)(int,int,double*,double*,double*,double*);
GetDP_Begin("Cal_deRhamTermOfFemEquation");
FI = EquationTerm_P->Case.LocalTerm.Active ;
QuantityStorageEqu_P = FI->QuantityStorageEqu_P ;
QuantityStorageDof_P = FI->QuantityStorageDof_P ;
/* ---------------------------------------------------------------------- */
/* E q u a t i o n s t o b u i l d ? */
/* ---------------------------------------------------------------------- */
/* Msg(INFO, "de Rham for Element %d", Element->Num); */
if (!(Nbr_Equ = QuantityStorageEqu_P->NbrElementaryBasisFunction)) {
if(Flag_VERBOSE > 2)
Msg(WARNING, "No equation to build in Element %d (de Rham)", Element->Num);
GetDP_End ;
}
/* ---------------------------------------------------------------------- */
/* 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 */
/* ---------------------------------------------------------------------- */
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 D e R h a m S u p p o r t */
/* ------------------------------------------------------------------- */
Get_NodesCoordinatesOfElement(Element) ;
Group_P = (struct Group*) List_Pointer(Problem_S.Group,
EquationTerm_P->Case.LocalTerm.InIndex) ;
Get_deRhamCells(Element, QuantityStorageEqu_P, Group_P, &Nbr_Cells,
Cells, Cell_Vector, &Cell_RelativeJacobianType);
/* All the cells created from one element are supposed to be of the same type */
/* ------------------------------------------------------------------- */
/* G e t J a c o b i a n M e t h o d f o r E l e m e n t */
/* ------------------------------------------------------------------- */
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) ;
/* ------------------------------------------------------------------- */
/* G e t J a c o b i a n M e t h o d f o r C e l l s */
/* ------------------------------------------------------------------- */
if(Cells[0].Type != POINT){
Cell_Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
Get_JacobianFunction(Element->JacobianCase->TypeJacobian +
Cell_RelativeJacobianType,
Cells[0].Type, &Cell_Type_Dimension) ;
}
/* ------------------------------------------------------------------------ */
/* C o m p u t a t i o n o f E l e m e n t a r y S u b m a t r i x */
/* ------------------------------------------------------------------------ */
/* Loop on source elements (> 1 only if integral quantity) */
while (1) {
if (FI->Type_DefineQuantityDof == INTEGRALQUANTITY) {
if ( Get_NextElementSource(Element->ElementSource) ) {
/* 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 ;
}
}
/* Elementary Submatrix initialization */
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. ;
/* Cell Value initialization */
for(i_Cell = 0 ; i_Cell < Nbr_Cells ; i_Cell ++)
for (j = 0 ; j < Nbr_Dof ; j++)
Cal_ZeroValue(&Cell_Value[i_Cell][j]);
/*
printf("Nb Equ= %d nbdof=%d cells =%d", Nbr_Equ , Nbr_Dof, Nbr_Cells);
*/
/* ------------------------------------- */
/* P o i n t c o l l o c a t i o n */
/* ------------------------------------- */
if(Cells[0].Type == POINT){
for(i_Cell = 0 ; i_Cell < Nbr_Cells ; i_Cell ++) {
Current.Element = Element ;
Current.u = Cells[i_Cell].x[0] ;
Current.v = Cells[i_Cell].y[0] ;
Current.w = Cells[i_Cell].z[0] ;
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) ;
/* Shape Functions (+ surrounding expression) */
Cal_vBFxDof(EquationTerm_P, FI,
QuantityStorage_P0, QuantityStorageDof_P,
Nbr_Dof, xFunctionBFDof, NULL, Cell_Value[i_Cell]);
}
}
/* ------------------------------------- */
/* I n t e g r a t i o n */
/* ------------------------------------- */
else {
switch (IntegrationCase_P->Type) {
case GAUSS :
case GAUSSLEGENDRE :
Quadrature_P = (struct Quadrature*)
List_PQuery(IntegrationCase_P->Case, &Cells[0].Type, fcmp_int);
if(!Quadrature_P)
Msg(GERROR, "Unknown type of Element (%s) for Integration method (%s)",
Get_StringForDefine(Element_Type, Cells[0].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_Cell = 0 ; i_Cell < Nbr_Cells ; i_Cell ++) {
for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
/* u,v,w in the Cell */
Get_IntPoint(Nbr_IntPoints, i_IntPoint, &u, &v, &w, &weight) ;
Get_BFGeoElement(&Cells[i_Cell], u, v, w) ;
Cells[i_Cell].DetJac = Cell_Get_Jacobian(&Cells[i_Cell], &Cells[i_Cell].Jac) ;
/* back to u,v,w in the original element. This is quite
inefficient (approximatively 3 jacobian inversions
instead of 1 if the transformation was defined in the
reverse order... -> to change ! */
x2 = y2 = z2 = 0. ;
for (i = 0 ; i < Cells[i_Cell].GeoElement->NbrNodes ; i++) {
x2 += Cells[i_Cell].x[i] * Cells[i_Cell].n[i];
y2 += Cells[i_Cell].y[i] * Cells[i_Cell].n[i];
z2 += Cells[i_Cell].z[i] * Cells[i_Cell].n[i];
}
xyz2uvwInAnElement(Element, x2, y2, z2, &u2, &v2, &w2) ;
Current.Element = Element ;
Current.u = u2 ;
Current.v = v2 ;
Current.w = w2 ;
/*
Msg(INFO, "Elm=%d uvw=%g %g %g xyz=%g %g %g",
Current.Element->Num, Current.u, Current.v, Current.w, x2, y2, z2);
*/
/* Shape Functions (+ surrounding expression) */
Cal_vBFxDof(EquationTerm_P, FI,
QuantityStorage_P0, QuantityStorageDof_P,
Nbr_Dof, xFunctionBFDof, NULL, vBFxDof);
Factor = (FI->Flag_ChangeCoord) ? weight * fabs(Cells[i_Cell].DetJac) : weight ;
/* Scalar product (by the tangent or normal to the cell). This is
definitely a hack and should be completely changed */
switch(FI->Type_FormDof){
case FORM0 :
case FORM3 :
case SCALAR :
break ;
case FORM1 :
case FORM2 :
case VECTOR :
for(j = 0 ; j < Nbr_Dof ; j++)
Cal_ProductValue(&Cell_Vector[i_Cell], &vBFxDof[j], &vBFxDof[j]) ;
break;
default :
Msg(GERROR, "Unknown type of Dof (%s) for scalar product in deRham Map",
Get_StringForDefine(Field_Type, FI->Type_FormDof));
break;
}
for (j = 0 ; j < Nbr_Dof ; j++)
for (k = 0 ; k < Current.NbrHar ; k++)
Cell_Value[i_Cell][j].Val[MAX_DIM*k] += Factor * vBFxDof[j].Val[MAX_DIM*k];
} /* for i_IntPoint ... */
} /* for i_Cell */
break ; /* case GAUSS */
default :
Msg(GERROR, "Unknown type of Integration method (%s)",
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
break;
} /* switch integration method */
} /* if point, else */
/* Discrete derivative if any */
switch (Group_P->FunctionType) {
case BOUNDARYOFDUALNODESOF :
Product = 1 ;
Dk = Geo_GetIM_Den_Xp(Element->Type, &Nbr_Ent1, &Nbr_Ent2); break;
break ;
case BOUNDARYOFDUALEDGESOF :
Product = 1 ;
Dk = Geo_GetIM_Dfe_Xp(Element->Type, &Nbr_Ent1, &Nbr_Ent2); break;
break;
default :
Product = 0 ;
break ;
}
/*
if(Nbr_Cells != Nbr_Ent1 || Nbr_Equ != Nbr_Ent2)
Msg(GERROR, "Error in deRham Map: wrong discrete derivative");
*/
if(Flag_VERBOSE == 10) {
printf("H D = ") ;
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("% .5e ", Cell_Value[i][j].Val[0]) ;
printf("]\n") ;
}
}
/* D^T (H D) if product ; H otherwise */
if(Product){
for (i = 0 ; i < Nbr_Equ ; i++)
for (j = 0 ; j < Nbr_Dof ; j++)
for (n = 0 ; n < Nbr_Cells ; n++)
if((ii = Dk [ n*Nbr_Ent1 + QuantityStorageEqu_P->BasisFunction[i].NumEntityInElement ])){
for (k = 0 ; k < Current.NbrHar ; k++)
Ek[i][j][k] += ii * Cell_Value[n][j].Val[MAX_DIM*k] ;
}
if(Flag_VERBOSE == 10) {
printf("D^T H D = ") ;
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("% .5e ", Ek[i][j][0]) ;
printf("]\n") ;
}
}
}
else{ /* ok, this is stupid */
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] = Cell_Value[i][j].Val[MAX_DIM*k] ;
}
/* 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