#define RCSID "$Id: Treatment_Formulation.c,v 1.19 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):
* David Colignon
* Johan Gyselinck
* Ruth Sabariego
*/
#include "GetDP.h"
#include "Treatment_Formulation.h"
#include "Get_DofOfElement.h"
#include "ExtendedGroup.h"
#include "GeoData.h"
#include "DofData.h"
#include "CurrentData.h"
#include "Tools.h"
#include "F_FMM.h"
#include "F_FMM_DTA.h"
#include "Pre_GetDofFMMGroup.h"
extern List_T * GeoData_L ;
int LastElementforAggreg;
/* ------------------------------------------------------------------------ */
/* T r e a t m e n t _ F o r m u l a t i o n */
/* ------------------------------------------------------------------------ */
void Treatment_Formulation(struct Formulation * Formulation_P) {
GetDP_Begin("Treatment_Formulation");
switch (Formulation_P->Type) {
case FEMEQUATION :
Treatment_FemFormulation(Formulation_P) ;
break ;
case GLOBALEQUATION :
Treatment_GlobalFormulation(Formulation_P) ;
break ;
default :
Msg(GERROR, "Unknown type for Formulation '%s'", Formulation_P->Name) ;
break ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* T r e a t m e n t _ F e m F o r m u l a t i o n */
/* ------------------------------------------------------------------------ */
void Treatment_FemFormulation(struct Formulation * Formulation_P) {
struct Element Element ;
struct QuantityStorage * QuantityStorage_P0, * QuantityStorage_P ;
struct QuantityStorage QuantityStorage_S ;
struct Dof DofForNoDof_P [NBR_MAX_HARMONIC] ;
struct EquationTerm * EquationTerm_P0 , * EquationTerm_P ;
struct GlobalQuantity * GlobalQuantity_P ;
int Nbr_DefineQuantity ;
struct DefineQuantity * DefineQuantity_P0 , * DefineQuantity_P ;
List_T * FemLocalTermActive_L ;
struct FemLocalTermActive* FemLocalTermActive_S ;
List_T * QuantityStorage_L ;
int i, j, Nbr_Element, i_Element, Nbr_EquationTerm, i_EquTerm ;
int Index_DefineQuantity, TraceGroupIndex_DefineQuantity ;
int Nbr_UnusedEquation ;
List_T * InitialListInIndex_L ;
int Nbr_Region, i_Region, Num_Region ;
char FileFMM[MAX_FILE_NAME_LENGTH];
extern struct Group * Generate_Group ;
extern double ** MH_Moving_Matrix ;
gMatrix A ;
gVector b ;
int Flag_Only ;
GetDP_Begin("Treatment_FemFormulation");
/* --------------------------------------------------------------- */
/* 0. Initialization of an active zone (QuantityStorage) for each */
/* DefineQuantity */
/* --------------------------------------------------------------- */
if (!(Nbr_EquationTerm = List_Nbr(Formulation_P->Equation)))
Msg(GERROR, "No equation in Formulation '%s'", Formulation_P->Name);
if (!(Nbr_DefineQuantity = List_Nbr(Formulation_P->DefineQuantity)))
Msg(GERROR, "No Quantity in Formulation '%s'", Formulation_P->Name);
DefineQuantity_P0 = (struct DefineQuantity*)
List_Pointer(Formulation_P->DefineQuantity, 0) ;
QuantityStorage_L = List_Create(Nbr_DefineQuantity, 1,
sizeof (struct QuantityStorage) ) ;
QuantityStorage_S.NumLastElementForFunctionSpace =
QuantityStorage_S.NumLastElementForDofDefinition =
QuantityStorage_S.NumLastElementForEquDefinition = 0 ;
for (i = 0 ; i < Nbr_DefineQuantity ; i++) {
QuantityStorage_S.DefineQuantity = DefineQuantity_P0 + i ;
if(QuantityStorage_S.DefineQuantity->Type == INTEGRALQUANTITY &&
QuantityStorage_S.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof < 0){
QuantityStorage_S.FunctionSpace = NULL ;
QuantityStorage_S.TypeQuantity = VECTOR ; /* to change */
}
else{
QuantityStorage_S.FunctionSpace = (struct FunctionSpace*)
List_Pointer(Problem_S.FunctionSpace,
QuantityStorage_S.DefineQuantity->FunctionSpaceIndex) ;
QuantityStorage_S.TypeQuantity = QuantityStorage_S.FunctionSpace->Type ;
}
List_Add(QuantityStorage_L, &QuantityStorage_S) ;
}
QuantityStorage_P0 = (struct QuantityStorage*)List_Pointer(QuantityStorage_L, 0) ;
Get_InitDofOfElement(&Element) ;
/* --------------------------------------------------------------- */
/* 1. Initialization of equation terms */
/* --------------------------------------------------------------- */
EquationTerm_P0 = (struct EquationTerm*)List_Pointer(Formulation_P->Equation, 0) ;
FemLocalTermActive_L = List_Create(Nbr_EquationTerm, 1,
sizeof (struct FemLocalTermActive) ) ;
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
List_Add(FemLocalTermActive_L, &FemLocalTermActive_S) ;
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
switch(EquationTerm_P->Type){
case GALERKIN :
EquationTerm_P->Case.LocalTerm.Active = (struct FemLocalTermActive*)
List_Pointer(FemLocalTermActive_L, i_EquTerm) ;
switch (TreatmentStatus) {
case _PRE :
Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
break ;
case _CAL :
Cal_InitGalerkinTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,
&QuantityStorage_S, DofForNoDof_P) ;
break ;
}
break;
case DERHAM :
EquationTerm_P->Case.LocalTerm.Active = (struct FemLocalTermActive*)
List_Pointer(FemLocalTermActive_L, i_EquTerm) ;
switch (TreatmentStatus) {
case _PRE :
Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
break ;
case _CAL :
Cal_InitdeRhamTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,
&QuantityStorage_S, DofForNoDof_P) ;
break ;
}
break;
case GLOBALTERM :
switch (TreatmentStatus) {
case _PRE :
Pre_InitGlobalTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;
break ;
}
break;
case GLOBALEQUATION :
break ;
default :
Msg(GERROR, "Unknown type of equation term") ;
break ;
}
}
/* ---------------------------------------------------------- */
/* PREprocessing FMM: Assigning NumDof to Groups */
/* --------------------------------------------------------- */
if(TreatmentStatus==_CAL && Flag_FMM == 1){
Msg(DIRECT, "P r e - P r o c e s s i n g F M M. . .") ;
Pre_GetDofFMMGroup(Formulation_P, QuantityStorage_P0) ;
if (!Current.FMM.NbrDir){ Msg(INFO, "No FMM far groups -> Flag_FMM = 0 ") ; Flag_FMM = 0 ; }
strcpy(FileFMM, Name_Generic); strcat(FileFMM, ".fmm"); Print_FMMGroupInfo(FileFMM);
Msg(DIRECT, "E n d P r e - P r o c e s s i n g F M M. . .") ;
}
/* ---------------------------------------------------------- */
/* 2. Loop on geometrical elements : */
/* Treatment of eventual GALERKIN / DERAHM terms */
/* --------------------------------------------------------- */
Nbr_Element = Geo_GetNbrGeoElements() ;
for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
if (Generate_Group) {
Element.Region = Geo_GetGeoElement(i_Element)->Region ;
while (i_Element < Nbr_Element &&
!List_Search(Generate_Group->InitialList,
&Element.Region, fcmp_int) ) {
i_Element++ ;
if (i_Element < Nbr_Element) Element.Region = Geo_GetGeoElement(i_Element)->Region ;
}
if (i_Element == Nbr_Element) break ;
}
Progress(i_Element, Nbr_Element, "") ;
Element.GeoElement = Geo_GetGeoElement(i_Element) ;
Element.Num = Element.GeoElement->Num ;
Element.Type = Element.GeoElement->Type ;
Current.Region = Element.Region = Element.GeoElement->Region ;
if (Flag_FMM) Element.FMMGroup = Element.GeoElement->FMMGroup ;
/* ---------------------------- */
/* 2.1. Loop on equation terms */
/* ---------------------------- */
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
if (Flag_FMM){
Current.FMM.Src = EquationTerm_P->Case.LocalTerm.FMMSource ;
Current.FMM.Obs = EquationTerm_P->Case.LocalTerm.FMMObservation ;
}
if (EquationTerm_P->Type == GALERKIN ||
EquationTerm_P->Type == DERHAM) {
/* if the element is in the support of integration of the term */
/* ----------------------------------------------------------- */
if (List_Search(((struct Group *)
List_Pointer(Problem_S.Group,
EquationTerm_P->Case.
LocalTerm.InIndex))->InitialList,
&Element.Region, fcmp_int ) ) {
if(Flag_VERBOSE == 10)
printf("==> Element #%d, EquationTerm #%d/%d\n",
Element.Num, i_EquTerm+1, Nbr_EquationTerm) ;
Current.IntegrationSupportIndex = EquationTerm_P->Case.LocalTerm.InIndex ;
/* ---------------------------------------------------------- */
/* 2.1.1. Loop on quantities (test fcts and shape functions) */
/* ---------------------------------------------------------- */
for (i = 0 ; i < EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex ; i++) {
Index_DefineQuantity =
EquationTerm_P->Case.LocalTerm.Term.QuantityIndexTable[i] ;
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
TraceGroupIndex_DefineQuantity =
EquationTerm_P->Case.LocalTerm.Term.QuantityTraceGroupIndexTable[i] ;
/* Only one analysis for each function space */
/*
* Attention : l'operateur de trace ne fonctionne que si le champ
* dont on prend la trace n'intervient qu'une seule fois dans le terme.
* du a - manque de generalite du code au niveau de la gestion des
* espaces fonctionnels des fcts tests pour 'Trace de Dof'
* - et Christophe fatigué
*/
if (QuantityStorage_P->NumLastElementForFunctionSpace != Element.Num ||
TraceGroupIndex_DefineQuantity >= 0) {
QuantityStorage_P->NumLastElementForFunctionSpace = Element.Num ;
switch (DefineQuantity_P->Type) {
case LOCALQUANTITY :
if(TraceGroupIndex_DefineQuantity >= 0){
Get_ElementTrace(&Element, TraceGroupIndex_DefineQuantity) ;
QuantityStorage_P->NumLastElementForFunctionSpace = Element.ElementTrace->Num ;
Get_DofOfElement
(Element.ElementTrace, QuantityStorage_P->FunctionSpace, QuantityStorage_P,
DefineQuantity_P->IndexInFunctionSpace) ;
}
else{
Get_DofOfElement
(&Element, QuantityStorage_P->FunctionSpace, QuantityStorage_P,
DefineQuantity_P->IndexInFunctionSpace) ;
}
break ;
case INTEGRALQUANTITY :
QuantityStorage_P->NbrElementaryBasisFunction = 0 ;
break ;
default :
Msg(GERROR, "Bad kind of Quantity in Formulation '%s'",
Formulation_P->Name);
break;
}
}
} /* for i = 0, 1 ... */
/* -------------------------------------- */
/* 2.1.2. Treatment of the equation term */
/* -------------------------------------- */
switch (TreatmentStatus) {
case _PRE :
Pre_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
break ;
case _PAR :
Par_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
break ;
case _CAL :
Flag_Only = 0 ;
if (Current.DofData->Flag_Only){
A = Current.DofData->A ;
b = Current.DofData->b ;
if (EquationTerm_P->Case.LocalTerm.MatrixIndex == -1)
EquationTerm_P->Case.LocalTerm.MatrixIndex = 0 ;
j = List_ISearch(Current.DofData->OnlyTheseMatrices,
&EquationTerm_P->Case.LocalTerm.MatrixIndex, fcmp_int) ;
if(j!=-1){
Flag_Only = 1 ;
switch(EquationTerm_P->Case.LocalTerm.MatrixIndex){
case 1 :
Current.DofData->A = Current.DofData->A1 ;
Current.DofData->b = Current.DofData->b1 ;
break;
case 2 :
Current.DofData->A = Current.DofData->A2 ;
Current.DofData->b = Current.DofData->b2 ;
break;
case 3 :
Current.DofData->A = Current.DofData->A3 ;
Current.DofData->b = Current.DofData->b3 ;
break;
}
}
}/* Only the matrices that vary are recalculated */
if (!Current.DofData->Flag_Only || (Current.DofData->Flag_Only && Flag_Only) ){
Nbr_UnusedEquation = 0;
QuantityStorage_P = QuantityStorage_P0 +
EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ;
if(Current.NbrCpu > 1){
for(i = 0 ; i < QuantityStorage_P->NbrElementaryBasisFunction ; i++){
if(QuantityStorage_P->BasisFunction[i].Dof->Type == DOF_UNKNOWN){
if(QuantityStorage_P->BasisFunction[i].Dof->Case.Unknown.NumDof >=
QuantityStorage_P->FunctionSpace->DofData->Part[Current.RankCpu] &&
QuantityStorage_P->BasisFunction[i].Dof->Case.Unknown.NumDof <
QuantityStorage_P->FunctionSpace->DofData->Part[Current.RankCpu+1])
break;
else
Nbr_UnusedEquation += 1;
}
}
}
if(Nbr_UnusedEquation != QuantityStorage_P->NbrElementaryBasisFunction){
if(EquationTerm_P->Type == GALERKIN)
Cal_GalerkinTermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
else
Cal_deRhamTermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
}
if (Current.DofData->Flag_Only && Flag_Only){
Current.DofData->A = A ;
Current.DofData->b = b ;
}
}/* Flag_Only */
break ;
case _CST :
Cst_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;
break ;
}
}/* if Support ... */
} /* if GALERKIN ... */
} /* for i_EquTerm ... */
} /* for i_Element ... */
if (MH_Moving_Matrix) {
List_Delete(FemLocalTermActive_L) ; List_Delete(QuantityStorage_L) ;
GetDP_End ;
}
/* ------------------------------------------------------ */
/* 3. Loop on equation terms : */
/* Treatment of eventual GLOBAL terms */
/* ------------------------------------------------------ */
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
if (EquationTerm_P->Type == GLOBALTERM) {
InitialListInIndex_L =
((struct Group *)List_Pointer(Problem_S.Group,
EquationTerm_P->Case.GlobalTerm.InIndex))
->InitialList ;
List_Sort(InitialListInIndex_L, fcmp_int) ;
Nbr_Region = List_Nbr(InitialListInIndex_L) ;
/* ---------------------------------------------- */
/* 3.1. Loop on Regions belonging to the support */
/* ---------------------------------------------- */
for (i_Region = 0 ; i_Region < Nbr_Region ; i_Region++) {
List_Read(InitialListInIndex_L, i_Region, &Num_Region) ;
Current.Region = Num_Region ;
/* ---------------------------------------------------------------- */
/* 3.1.1. Loop on Quantities (test functions and shape functions) */
/* ---------------------------------------------------------------- */
for (i = 0 ; i < EquationTerm_P->Case.GlobalTerm.Term.NbrQuantityIndex ; i++) {
Index_DefineQuantity =
EquationTerm_P->Case.GlobalTerm.Term.QuantityIndexTable[i] ;
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
GlobalQuantity_P = (struct GlobalQuantity*)
List_Pointer(QuantityStorage_P->FunctionSpace->GlobalQuantity,
*(int*)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
/* Only one Function space analysis */
/* -------------------------------- */
if (QuantityStorage_P->NumLastElementForFunctionSpace != Num_Region) {
QuantityStorage_P->NumLastElementForFunctionSpace = Num_Region ;
switch (DefineQuantity_P->Type) {
case GLOBALQUANTITY :
Get_DofOfRegion
(Num_Region, GlobalQuantity_P,
QuantityStorage_P->FunctionSpace, QuantityStorage_P) ;
break ;
default :
Msg(GERROR, "Bad kind of Quantity in Formulation '%s'",
Formulation_P->Name);
break;
}
}
} /* for i = 0, 1 ... */
/* ------------------------------ */
/* 3.1.2. Treatment of the term */
/* ------------------------------ */
switch (TreatmentStatus) {
case _PRE :
Pre_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
QuantityStorage_P0) ;
break ;
case _PAR :
Msg(GERROR, "Partitionning not done for Global Terms");
break;
case _CAL :
Cal_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
QuantityStorage_P0,
&QuantityStorage_S, DofForNoDof_P) ;
break ;
}
}
} /* if GLOBALTERM ... */
} /* for i_EquTerm ... */
/* --------------------------------------------------------- */
/* 4. Loop on equation terms : */
/* Treatment of eventual GLOBAL EQUATION terms */
/* --------------------------------------------------------- */
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
if (EquationTerm_P->Type == GLOBALEQUATION) {
if (EquationTerm_P->Case.GlobalEquation.ConstraintIndex >= 0)
switch (TreatmentStatus) {
case _PRE :
Pre_FemGlobalEquation(EquationTerm_P,
DefineQuantity_P0, QuantityStorage_P0) ;
break ;
case _PAR :
Msg(GERROR, "Partitionning not done for Global Equations");
break;
case _CAL :
Cal_FemGlobalEquation(EquationTerm_P,
DefineQuantity_P0, QuantityStorage_P0) ;
break ;
}
} /* if GLOBALEQUATION ... */
} /* for i_EquTerm ... */
if (TreatmentStatus==_CAL && Flag_FMM ){
Msg(DIRECT, "P r o c e s s i n g F M M D A T A. . .") ;
Cal_FMMGalerkinAggregation(EquationTerm_P0, QuantityStorage_P0) ;
Cal_FMMGalerkinDisaggregation(EquationTerm_P0, QuantityStorage_P0) ;
if (Flag_FMM != 1) Msg(DIRECT, "E n d P r o c e s s i n g F M M D A T A. . .") ;
}
if (TreatmentStatus==_CAL && Flag_FMM == 1){
if (Current.FMM.Precision == 0.)/* Helmholtz */
GF_FMMTranslationValue() ;
else
GF_FMMTranslationValueAdaptive() ;
Msg(DIRECT, "E n d P r o c e s s i n g F M M D A T A. . .") ;
Flag_FMM = 2 ;
} /* FMM */
/* -------------------------- */
/* 5. End of FEM treatment */
/* -------------------------- */
List_Delete(FemLocalTermActive_L) ; List_Delete(QuantityStorage_L) ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* C a l _ F e m G l o b a l E q u a t i o n */
/* ------------------------------------------------------------------------ */
void Cal_FemGlobalEquation(struct EquationTerm * EquationTerm_P,
struct DefineQuantity * DefineQuantity_P0,
struct QuantityStorage * QuantityStorage_P0) {
int Nbr_GlobalEquationTerm, i_GlobalEquationTerm ;
struct Constraint * Constraint_P ;
struct GlobalEquationTerm * GlobalEquationTerm_P ;
int Nbr_EquAndDof ;
List_T * InitialListInIndex_L, * RegionIndex_L ;
int Nbr_Region, i_Region, Num_Region, k ;
int Nbr_MCPR, i_MCPR, Nbr_CPR, i_CPR, i_Node, i_Loop, j_Branch, Num_Equ ;
struct MultiConstraintPerRegion * MCPR_P ;
struct ConstraintPerRegion * CPR_P ;
struct Group * Group_P ;
double Val[NBR_MAX_HARMONIC] ;
struct DofGlobal { int NumRegion ; struct Dof * Dof ; } ;
List_T * DofGlobal_Equ_L, * DofGlobal_DofNode_L, * DofGlobal_DofLoop_L ;
struct DofGlobal * DofGlobal_Equ, * DofGlobal_DofNode, * DofGlobal_DofLoop ;
struct DofGlobal DofGlobal_S, * DofGlobal_P ;
struct Dof * Cal_FemGlobalEquation2(int Index_DefineQuantity, int Num_Region,
struct DefineQuantity * DefineQuantity_P0,
struct QuantityStorage * QuantityStorage_P0) ;
GetDP_Begin("Cal_FemGlobalEquation");
/* Liste des Regions auxquelles on associe des Equations de Type 'Network' */
RegionIndex_L = List_Create(50,50, sizeof(int)) ;
Constraint_P = (struct Constraint*)
List_Pointer(Problem_S.Constraint,
EquationTerm_P->Case.GlobalEquation.ConstraintIndex) ;
Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion) ;
for (i_MCPR = 0 ; i_MCPR < Nbr_MCPR ; i_MCPR++) {
MCPR_P = (struct MultiConstraintPerRegion*)
List_Pointer(Constraint_P->MultiConstraintPerRegion, i_MCPR) ;
Nbr_CPR = List_Nbr(MCPR_P->ConstraintPerRegion) ;
for (i_CPR = 0 ; i_CPR < Nbr_CPR ; i_CPR++) {
CPR_P = (struct ConstraintPerRegion*)
List_Pointer(MCPR_P->ConstraintPerRegion, i_CPR) ;
Group_P = (struct Group *)
List_Pointer(Problem_S.Group, CPR_P->RegionIndex) ;
List_Read(Group_P->InitialList, 0, &Num_Region) ;
if (!List_Search(RegionIndex_L, &Num_Region, fcmp_int))
List_Add(RegionIndex_L, &Num_Region) ;
else {
Msg(GERROR, "2 occurences of Elementary Region #%d in Contraint '%s'",
Num_Region, Constraint_P->Name);
}
}
}
Nbr_EquAndDof = List_Nbr(RegionIndex_L) ;
if (!Nbr_EquAndDof){
GetDP_End ;
}
DofGlobal_Equ_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal)) ;
DofGlobal_DofNode_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal)) ;
DofGlobal_DofLoop_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal)) ;
/* Construction des listes de Dof globaux pour Equ, DofNode, DofLoop */
Nbr_GlobalEquationTerm =
List_Nbr(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm) ;
for (i_GlobalEquationTerm = 0 ;
i_GlobalEquationTerm < Nbr_GlobalEquationTerm ; i_GlobalEquationTerm++) {
GlobalEquationTerm_P = (struct GlobalEquationTerm*)
List_Pointer(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm,
i_GlobalEquationTerm) ;
InitialListInIndex_L =
((struct Group *)List_Pointer(Problem_S.Group,
GlobalEquationTerm_P->InIndex))->InitialList ;
Nbr_Region = List_Nbr(InitialListInIndex_L) ;
List_Sort(InitialListInIndex_L, fcmp_int) ;
for (i_Region = 0 ; i_Region < Nbr_Region ; i_Region++) {
List_Read(InitialListInIndex_L, i_Region, &Num_Region) ;
if (List_Search(RegionIndex_L, &Num_Region, fcmp_int)) {
DofGlobal_S.NumRegion = Num_Region ;
DofGlobal_S.Dof = Cal_FemGlobalEquation2
(GlobalEquationTerm_P->DefineQuantityIndexEqu, Num_Region,
DefineQuantity_P0, QuantityStorage_P0) ;
List_Add(DofGlobal_Equ_L, &DofGlobal_S) ;
DofGlobal_S.Dof = Cal_FemGlobalEquation2
(GlobalEquationTerm_P->DefineQuantityIndexNode, Num_Region,
DefineQuantity_P0, QuantityStorage_P0) ;
List_Add(DofGlobal_DofNode_L, &DofGlobal_S) ;
DofGlobal_S.Dof = Cal_FemGlobalEquation2
(GlobalEquationTerm_P->DefineQuantityIndexLoop, Num_Region,
DefineQuantity_P0, QuantityStorage_P0) ;
List_Add(DofGlobal_DofLoop_L, &DofGlobal_S) ;
}
}
}
if (List_Nbr(DofGlobal_Equ_L) != Nbr_EquAndDof) {
Msg(DIRECT, ERROR_STR "Incompatible number of equations with Contraint '%s'.",
Constraint_P->Name);
Msg(GERROR, "(%d equations obtained while %d branches are defined)",
List_Nbr(DofGlobal_Equ_L), Nbr_EquAndDof);
}
DofGlobal_Equ = (struct DofGlobal*)List_Pointer(DofGlobal_Equ_L , 0) ;
DofGlobal_DofNode = (struct DofGlobal*)List_Pointer(DofGlobal_DofNode_L, 0) ;
DofGlobal_DofLoop = (struct DofGlobal*)List_Pointer(DofGlobal_DofLoop_L, 0) ;
for (k = 0 ; k < List_Nbr(DofGlobal_Equ_L) ; k++) {
if (DofGlobal_Equ[k].Dof->Type == DOF_FIXED ||
DofGlobal_Equ[k].Dof->Type == DOF_LINK) {
if (DofGlobal_Equ[k].Dof == DofGlobal_DofNode[k].Dof)
DofGlobal_Equ[k].Dof = DofGlobal_DofLoop[k].Dof ;
else
DofGlobal_Equ[k].Dof = DofGlobal_DofNode[k].Dof ;
}
}
/* Construction des equations (assemblage) */
Num_Equ = 0 ;
Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion) ;
for (i_MCPR = 0 ; i_MCPR < Nbr_MCPR ; i_MCPR++) {
MCPR_P = (struct MultiConstraintPerRegion*)
List_Pointer(Constraint_P->MultiConstraintPerRegion, i_MCPR) ;
if (!MCPR_P->Active)
MCPR_P->Active = Generate_Network(MCPR_P->ConstraintPerRegion) ;
for (i_Node = 0 ; i_Node < MCPR_P->Active->Case.Network.NbrNode ; i_Node++) {
for (j_Branch = 0 ;
j_Branch < MCPR_P->Active->Case.Network.NbrBranch ; j_Branch++) {
if (MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]) {
Group_P = (struct Group *)
List_Pointer
(Problem_S.Group,
((struct ConstraintPerRegion *)
List_Pointer(MCPR_P->ConstraintPerRegion, j_Branch))->RegionIndex) ;
List_Read(Group_P->InitialList, 0, &Num_Region) ;
DofGlobal_P = (struct DofGlobal*)
List_PQuery(DofGlobal_DofNode_L, &Num_Region, fcmp_int) ;
Val[0] =
(double)(MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]) ;
if (Current.NbrHar > 1) {
Val[1] = 0. ;
for (k = 2 ; k < Current.NbrHar ; k += 2) {
Val[k] = Val[0] ; Val[k+1] = 0. ;
}
}
/*
Msg(CHECK, "Node: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",
Num_Equ,
DofGlobal_Equ[Num_Equ].Dof->NumType, DofGlobal_Equ[Num_Equ].Dof->Entity,
DofGlobal_P->Dof->NumType, DofGlobal_P->Dof->Entity,
Val[0]
) ;
*/
Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof, Val) ;
}
}
Num_Equ++ ;
} /* for i_Node ... */
for (i_Loop = 0 ; i_Loop < MCPR_P->Active->Case.Network.NbrLoop ; i_Loop++) {
for (j_Branch = 0 ;
j_Branch < MCPR_P->Active->Case.Network.NbrBranch ; j_Branch++) {
if (MCPR_P->Active->Case.Network.MatLoop[i_Loop][j_Branch]) {
Group_P = (struct Group *)
List_Pointer
(Problem_S.Group,
((struct ConstraintPerRegion *)
List_Pointer(MCPR_P->ConstraintPerRegion, j_Branch))->RegionIndex) ;
List_Read(Group_P->InitialList, 0, &Num_Region) ;
DofGlobal_P = (struct DofGlobal*)
List_PQuery(DofGlobal_DofLoop_L, &Num_Region, fcmp_int) ;
Val[0] =
(double)(MCPR_P->Active->Case.Network.MatLoop[i_Loop][j_Branch]) ;
if (Current.NbrHar > 1) {
Val[1] = 0. ;
for (k = 2 ; k < Current.NbrHar ; k += 2) {
Val[k] = Val[0] ; Val[k+1] = 0. ;
}
}
/*
Msg(CHECK, "Loop: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",
Num_Equ,
DofGlobal_Equ[Num_Equ].Dof->NumType, DofGlobal_Equ[Num_Equ].Dof->Entity,
DofGlobal_P->Dof->NumType, DofGlobal_P->Dof->Entity,
Val[0]
) ;
*/
Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof, Val) ;
}
}
Num_Equ++ ;
} /* for i_Loop ... */
} /* for i_MCPR ... */
List_Delete(DofGlobal_Equ_L) ;
List_Delete(DofGlobal_DofNode_L) ; List_Delete(DofGlobal_DofLoop_L) ;
List_Delete(RegionIndex_L) ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
struct Dof * Cal_FemGlobalEquation2(int Index_DefineQuantity, int Num_Region,
struct DefineQuantity * DefineQuantity_P0,
struct QuantityStorage * QuantityStorage_P0) {
struct DefineQuantity * DefineQuantity_P ;
struct QuantityStorage * QuantityStorage_P ;
struct GlobalQuantity * GlobalQuantity_P ;
struct QuantityStorage QuaSto_S ;
GetDP_Begin("Cal_FemGlobalEquation2");
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;
GlobalQuantity_P = (struct GlobalQuantity*)
List_Pointer(QuantityStorage_P->FunctionSpace->GlobalQuantity,
*(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0)) ;
Get_DofOfRegion(Num_Region, GlobalQuantity_P,
QuantityStorage_P->FunctionSpace, &QuaSto_S) ;
if (QuaSto_S.NbrElementaryBasisFunction == 1){
GetDP_Return(QuaSto_S.BasisFunction[0].Dof) ;
}
else {
Msg(GERROR, "Not 1 Dof associated with GlobalQuantity (Region #%d)", Num_Region) ;
GetDP_Return(NULL) ;
}
}
/* ------------------------------------------------------------------------ */
/* T r e a t m e n t _ G l o b a l F o r m u l a t i o n */
/* ------------------------------------------------------------------------ */
void Treatment_GlobalFormulation(struct Formulation * Formulation_P) {
GetDP_Begin("Treatment_GlobalFormulation");
Msg(GERROR, "You should not be here!") ;
GetDP_End ;
}
syntax highlighted by Code2HTML, v. 0.9.1