"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Kernel/Treatment_Formulation.cpp" between
getdp-3.4.0-source.tgz and getdp-3.5.0-source.tgz

About: GetDP is a general finite element solver using mixed elements to discretize de Rham-type complexes in one, two and three dimensions.

Treatment_Formulation.cpp  (getdp-3.4.0-source.tgz)	:	Treatment_Formulation.cpp  (getdp-3.5.0-source.tgz)
skipping to change at line 20		skipping to change at line 20
#include "Get_ElementSource.h"		#include "Get_ElementSource.h"
#include "Pre_TermOfFemEquation.h"		#include "Pre_TermOfFemEquation.h"
#include "Cal_GalerkinTermOfFemEquation.h"		#include "Cal_GalerkinTermOfFemEquation.h"
#include "Cal_SmallFemTermOfFemEquation.h"		#include "Cal_SmallFemTermOfFemEquation.h"
#include "Cal_GlobalTermOfFemEquation.h"		#include "Cal_GlobalTermOfFemEquation.h"
#include "Cal_AssembleTerm.h"		#include "Cal_AssembleTerm.h"
#include "Generate_Network.h"		#include "Generate_Network.h"
#include "ExtendedGroup.h"		#include "ExtendedGroup.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;
extern int TreatmentStatus ;		extern int TreatmentStatus;
extern List_T * GeoData_L ;		extern List_T *GeoData_L;
/* ------------------------------------------------------------------------ */		/* ------------------------------------------------------------------------ */
/* C a l _ F e m G l o b a l E q u a t i o n */		/* C a l _ F e m G l o b a l E q u a t i o n */
/* ------------------------------------------------------------------------ */		/* ------------------------------------------------------------------------ */
struct Dof * Cal_FemGlobalEquation2(int Index_DefineQuantity, int Num_Region,		struct Dof *Cal_FemGlobalEquation2(int Index_DefineQuantity, int Num_Region,
struct DefineQuantity * DefineQuantity_P0,		struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage * QuantityStorage_P0)		struct QuantityStorage *QuantityStorage_P0)
{		{
struct DefineQuantity * DefineQuantity_P ;		struct DefineQuantity *DefineQuantity_P;
struct QuantityStorage * QuantityStorage_P ;		struct QuantityStorage *QuantityStorage_P;
struct GlobalQuantity * GlobalQuantity_P ;		struct GlobalQuantity *GlobalQuantity_P;
struct QuantityStorage QuaSto_S ;		struct QuantityStorage QuaSto_S;
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;		DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;		QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity;
GlobalQuantity_P = (struct GlobalQuantity*)		GlobalQuantity_P = (struct GlobalQuantity *)List_Pointer(
List_Pointer(QuantityStorage_P->FunctionSpace->GlobalQuantity,		QuantityStorage_P->FunctionSpace->GlobalQuantity,
*(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0))		*(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace, 0));
;
Get_DofOfRegion(Num_Region, GlobalQuantity_P,		Get_DofOfRegion(Num_Region, GlobalQuantity_P,
QuantityStorage_P->FunctionSpace, &QuaSto_S) ;		QuantityStorage_P->FunctionSpace, &QuaSto_S);
if (QuaSto_S.NbrElementaryBasisFunction == 1){		if(QuaSto_S.NbrElementaryBasisFunction == 1) {
return QuaSto_S.BasisFunction[0].Dof ;		return QuaSto_S.BasisFunction[0].Dof;
}		}
else {		else {
Message::Error( "Not 1 Dof associated with GlobalQuantity (Region #%d)", Num		Message::Error("Not 1 Dof associated with GlobalQuantity (Region #%d)",
_Region) ;		Num_Region);
return NULL ;		return NULL;
}		}
}		}
void Cal_FemGlobalEquation(struct EquationTerm * EquationTerm_P,		void Cal_FemGlobalEquation(struct EquationTerm *EquationTerm_P,
struct DefineQuantity * DefineQuantity_P0,		struct DefineQuantity *DefineQuantity_P0,
struct QuantityStorage * QuantityStorage_P0)		struct QuantityStorage *QuantityStorage_P0)
{		{
int Nbr_GlobalEquationTerm, i_GlobalEquationTerm ;		int Nbr_GlobalEquationTerm, i_GlobalEquationTerm;
struct Constraint * Constraint_P ;		struct Constraint *Constraint_P;
struct GlobalEquationTerm * GlobalEquationTerm_P ;		struct GlobalEquationTerm *GlobalEquationTerm_P;
int Nbr_EquAndDof ;		int Nbr_EquAndDof;
List_T * InitialListInIndex_L, * RegionIndex_L ;		List_T *InitialListInIndex_L, *RegionIndex_L;
int Nbr_Region, i_Region, Num_Region, k ;		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 ;		int Nbr_MCPR, i_MCPR, Nbr_CPR, i_CPR, i_Node, i_Loop, j_Branch, Num_Equ;
struct MultiConstraintPerRegion * MCPR_P ;		struct MultiConstraintPerRegion *MCPR_P;
struct ConstraintPerRegion * CPR_P ;		struct ConstraintPerRegion *CPR_P;
struct Group * Group_P ;		struct Group *Group_P;
double Val[NBR_MAX_HARMONIC] ;		double Val[NBR_MAX_HARMONIC];
struct DofGlobal { int NumRegion ; struct Dof * Dof ; } ;		struct DofGlobal {
List_T * DofGlobal_Equ_L, * DofGlobal_DofNode_L, * DofGlobal_DofLoop_L ;		int NumRegion;
struct DofGlobal * DofGlobal_Equ, * DofGlobal_DofNode, * DofGlobal_DofLoop ;		struct Dof *Dof;
struct DofGlobal DofGlobal_S, * DofGlobal_P ;		};
		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;
/* Liste des Regions auxquelles on associe des Equations de Type 'Network' */		/* Liste des Regions auxquelles on associe des Equations de Type 'Network' */
RegionIndex_L = List_Create(50,50, sizeof(int)) ;		RegionIndex_L = List_Create(50, 50, sizeof(int));
Constraint_P = (struct Constraint*)		Constraint_P = (struct Constraint *)List_Pointer(
List_Pointer(Problem_S.Constraint,		Problem_S.Constraint, EquationTerm_P->Case.GlobalEquation.ConstraintIndex);
EquationTerm_P->Case.GlobalEquation.ConstraintIndex) ;		Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion);
Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion) ;		for(i_MCPR = 0; i_MCPR < Nbr_MCPR; i_MCPR++) {
for (i_MCPR = 0 ; i_MCPR < Nbr_MCPR ; i_MCPR++) {		MCPR_P = (struct MultiConstraintPerRegion *)List_Pointer(
MCPR_P = (struct MultiConstraintPerRegion*)		Constraint_P->MultiConstraintPerRegion, i_MCPR);
List_Pointer(Constraint_P->MultiConstraintPerRegion, i_MCPR) ;		Nbr_CPR = List_Nbr(MCPR_P->ConstraintPerRegion);
Nbr_CPR = List_Nbr(MCPR_P->ConstraintPerRegion) ;		for(i_CPR = 0; i_CPR < Nbr_CPR; i_CPR++) {
for (i_CPR = 0 ; i_CPR < Nbr_CPR ; i_CPR++) {		CPR_P = (struct ConstraintPerRegion *)List_Pointer(
CPR_P = (struct ConstraintPerRegion*)		MCPR_P->ConstraintPerRegion, i_CPR);
List_Pointer(MCPR_P->ConstraintPerRegion, i_CPR) ;		Group_P =
Group_P = (struct Group *)		(struct Group *)List_Pointer(Problem_S.Group, CPR_P->RegionIndex);
List_Pointer(Problem_S.Group, CPR_P->RegionIndex) ;		List_Read(Group_P->InitialList, 0, &Num_Region);
List_Read(Group_P->InitialList, 0, &Num_Region) ;		if(!List_Search(RegionIndex_L, &Num_Region, fcmp_int))
if (!List_Search(RegionIndex_L, &Num_Region, fcmp_int))		List_Add(RegionIndex_L, &Num_Region);
List_Add(RegionIndex_L, &Num_Region) ;
else {		else {
Message::Error("2 occurences of Elementary Region #%d in Contraint '%s'",		Message::Error(
Num_Region, Constraint_P->Name);		"2 occurences of Elementary Region #%d in Contraint '%s'", Num_Region,
		Constraint_P->Name);
return;		return;
}		}
}		}
}		}
Nbr_EquAndDof = List_Nbr(RegionIndex_L) ;		Nbr_EquAndDof = List_Nbr(RegionIndex_L);
if (!Nbr_EquAndDof){		if(!Nbr_EquAndDof) { return; }
return ;
}
DofGlobal_Equ_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal))		DofGlobal_Equ_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal));
;		DofGlobal_DofNode_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));
;
DofGlobal_DofLoop_L = List_Create(Nbr_EquAndDof, 1, sizeof(struct DofGlobal))
;
/* Construction des listes de Dof globaux pour Equ, DofNode, DofLoop */		/* Construction des listes de Dof globaux pour Equ, DofNode, DofLoop */
Nbr_GlobalEquationTerm =		Nbr_GlobalEquationTerm =
List_Nbr(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm) ;		List_Nbr(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm);
for (i_GlobalEquationTerm = 0 ;		for(i_GlobalEquationTerm = 0; i_GlobalEquationTerm < Nbr_GlobalEquationTerm;
i_GlobalEquationTerm < Nbr_GlobalEquationTerm ; i_GlobalEquationTerm++) {		i_GlobalEquationTerm++) {
GlobalEquationTerm_P = (struct GlobalEquationTerm*)		GlobalEquationTerm_P = (struct GlobalEquationTerm *)List_Pointer(
List_Pointer(EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm,		EquationTerm_P->Case.GlobalEquation.GlobalEquationTerm,
i_GlobalEquationTerm) ;		i_GlobalEquationTerm);
InitialListInIndex_L =		InitialListInIndex_L = ((struct Group *)List_Pointer(
((struct Group *)List_Pointer(Problem_S.Group,		Problem_S.Group, GlobalEquationTerm_P->InIndex))
GlobalEquationTerm_P->InIndex))->InitialList		->InitialList;
;		Nbr_Region = List_Nbr(InitialListInIndex_L);
Nbr_Region = List_Nbr(InitialListInIndex_L) ;		List_Sort(InitialListInIndex_L, fcmp_int);
List_Sort(InitialListInIndex_L, fcmp_int) ;
		for(i_Region = 0; i_Region < Nbr_Region; i_Region++) {
for (i_Region = 0 ; i_Region < Nbr_Region ; i_Region++) {		List_Read(InitialListInIndex_L, i_Region, &Num_Region);
List_Read(InitialListInIndex_L, i_Region, &Num_Region) ;		if(List_Search(RegionIndex_L, &Num_Region, fcmp_int)) {
if (List_Search(RegionIndex_L, &Num_Region, fcmp_int)) {		DofGlobal_S.NumRegion = Num_Region;
DofGlobal_S.NumRegion = Num_Region ;		DofGlobal_S.Dof = Cal_FemGlobalEquation2(
DofGlobal_S.Dof = Cal_FemGlobalEquation2		GlobalEquationTerm_P->DefineQuantityIndexEqu, Num_Region,
(GlobalEquationTerm_P->DefineQuantityIndexEqu, Num_Region,		DefineQuantity_P0, QuantityStorage_P0);
DefineQuantity_P0, QuantityStorage_P0) ;		List_Add(DofGlobal_Equ_L, &DofGlobal_S);
List_Add(DofGlobal_Equ_L, &DofGlobal_S) ;		DofGlobal_S.Dof = Cal_FemGlobalEquation2(
DofGlobal_S.Dof = Cal_FemGlobalEquation2		GlobalEquationTerm_P->DefineQuantityIndexNode, Num_Region,
(GlobalEquationTerm_P->DefineQuantityIndexNode, Num_Region,		DefineQuantity_P0, QuantityStorage_P0);
DefineQuantity_P0, QuantityStorage_P0) ;		List_Add(DofGlobal_DofNode_L, &DofGlobal_S);
List_Add(DofGlobal_DofNode_L, &DofGlobal_S) ;		DofGlobal_S.Dof = Cal_FemGlobalEquation2(
DofGlobal_S.Dof = Cal_FemGlobalEquation2		GlobalEquationTerm_P->DefineQuantityIndexLoop, Num_Region,
(GlobalEquationTerm_P->DefineQuantityIndexLoop, Num_Region,		DefineQuantity_P0, QuantityStorage_P0);
DefineQuantity_P0, QuantityStorage_P0) ;		List_Add(DofGlobal_DofLoop_L, &DofGlobal_S);
List_Add(DofGlobal_DofLoop_L, &DofGlobal_S) ;
}		}
}		}
}		}
if (List_Nbr(DofGlobal_Equ_L) != Nbr_EquAndDof) {		if(List_Nbr(DofGlobal_Equ_L) != Nbr_EquAndDof) {
Message::Error("Incompatible number of equations with Contraint '%s' "		Message::Error("Incompatible number of equations with Contraint '%s' "
"(%d equations obtained while %d branches are defined)",		"(%d equations obtained while %d branches are defined)",
Constraint_P->Name, List_Nbr(DofGlobal_Equ_L), Nbr_EquAndDof)		Constraint_P->Name, List_Nbr(DofGlobal_Equ_L),
;		Nbr_EquAndDof);
return;		return;
}		}
DofGlobal_Equ = (struct DofGlobal*)List_Pointer(DofGlobal_Equ_L , 0) ;		DofGlobal_Equ = (struct DofGlobal *)List_Pointer(DofGlobal_Equ_L, 0);
DofGlobal_DofNode = (struct DofGlobal*)List_Pointer(DofGlobal_DofNode_L, 0) ;		DofGlobal_DofNode = (struct DofGlobal *)List_Pointer(DofGlobal_DofNode_L, 0);
DofGlobal_DofLoop = (struct DofGlobal*)List_Pointer(DofGlobal_DofLoop_L, 0) ;		DofGlobal_DofLoop = (struct DofGlobal *)List_Pointer(DofGlobal_DofLoop_L, 0);
for (k = 0 ; k < List_Nbr(DofGlobal_Equ_L) ; k++) {		for(k = 0; k < List_Nbr(DofGlobal_Equ_L); k++) {
if (DofGlobal_Equ[k].Dof->Type == DOF_FIXED ||		if(DofGlobal_Equ[k].Dof->Type == DOF_FIXED ||
DofGlobal_Equ[k].Dof->Type == DOF_LINK) {		DofGlobal_Equ[k].Dof->Type == DOF_LINK) {
if (DofGlobal_Equ[k].Dof == DofGlobal_DofNode[k].Dof)		if(DofGlobal_Equ[k].Dof == DofGlobal_DofNode[k].Dof)
DofGlobal_Equ[k].Dof = DofGlobal_DofLoop[k].Dof ;		DofGlobal_Equ[k].Dof = DofGlobal_DofLoop[k].Dof;
else		else
DofGlobal_Equ[k].Dof = DofGlobal_DofNode[k].Dof ;		DofGlobal_Equ[k].Dof = DofGlobal_DofNode[k].Dof;
}		}
}		}
/* Construction des equations (assemblage) */		/* Construction des equations (assemblage) */
Num_Equ = 0 ;		Num_Equ = 0;
Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion) ;		Nbr_MCPR = List_Nbr(Constraint_P->MultiConstraintPerRegion);
for (i_MCPR = 0 ; i_MCPR < Nbr_MCPR ; i_MCPR++) {		for(i_MCPR = 0; i_MCPR < Nbr_MCPR; i_MCPR++) {
MCPR_P = (struct MultiConstraintPerRegion*)		MCPR_P = (struct MultiConstraintPerRegion *)List_Pointer(
List_Pointer(Constraint_P->MultiConstraintPerRegion, i_MCPR) ;		Constraint_P->MultiConstraintPerRegion, i_MCPR);
if (!MCPR_P->Active)		if(!MCPR_P->Active)
MCPR_P->Active = Generate_Network(MCPR_P->Name, MCPR_P->ConstraintPerRegio		MCPR_P->Active =
n) ;		Generate_Network(MCPR_P->Name, MCPR_P->ConstraintPerRegion);
for (i_Node = 0 ; i_Node < MCPR_P->Active->Case.Network.NbrNode ; i_Node++)		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;
for (j_Branch = 0 ;		j_Branch++) {
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(
if (MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]) {		Problem_S.Group, ((struct ConstraintPerRegion *)List_Pointer(
		MCPR_P->ConstraintPerRegion, j_Branch))
Group_P = (struct Group *)		->RegionIndex);
List_Pointer		List_Read(Group_P->InitialList, 0, &Num_Region);
(Problem_S.Group,
((struct ConstraintPerRegion *)		DofGlobal_P = (struct DofGlobal *)List_PQuery(DofGlobal_DofNode_L,
List_Pointer(MCPR_P->ConstraintPerRegion, j_Branch))->RegionIndex		&Num_Region, fcmp_int);
) ;
List_Read(Group_P->InitialList, 0, &Num_Region) ;		Val[0] =
		(double)(MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]);
DofGlobal_P = (struct DofGlobal*)		if(Current.NbrHar > 1) {
List_PQuery(DofGlobal_DofNode_L, &Num_Region, fcmp_int) ;		Val[1] = 0.;
		for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) {
Val[0] =		Val[k] = Val[0];
(double)(MCPR_P->Active->Case.Network.MatNode[i_Node][j_Branch]) ;		Val[k + 1] = 0.;
if (Current.NbrHar > 1) {		}
Val[1] = 0. ;		}
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2		/*
) {		Message::Info("Node: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",
Val[k] = Val[0] ; Val[k+1] = 0. ;		Num_Equ,
}		DofGlobal_Equ[Num_Equ].Dof->NumType,
}		DofGlobal_Equ[Num_Equ].Dof->Entity, DofGlobal_P->Dof->NumType,
/*		DofGlobal_P->Dof->Entity, Val[0] ) ;
Message::Info("Node: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",		*/
Num_Equ,		Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof,
DofGlobal_Equ[Num_Equ].Dof->NumType, DofGlobal_Equ[Num_Equ].Dof->En		Val);
tity,		}
DofGlobal_P->Dof->NumType, DofGlobal_P->Dof->Entity,
Val[0]
) ;
*/
Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof,
Val) ;
}
}		}
Num_Equ++ ;		Num_Equ++;
} /* for i_Node ... */		} /* 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 < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2
) {
Val[k] = Val[0] ; Val[k+1] = 0. ;
}
}
/*
Message::Info("Loop: eq.(%d) [%d, %d], dof [%d, %d] : %.16g\n",
Num_Equ,
DofGlobal_Equ[Num_Equ].Dof->NumType, DofGlobal_Equ[Num_Equ].Dof->En
tity,
DofGlobal_P->Dof->NumType, DofGlobal_P->Dof->Entity,
Val[0]
) ;
*/
Cal_AssembleTerm_NeverDt(DofGlobal_Equ[Num_Equ].Dof, DofGlobal_P->Dof,
Val) ;
}		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 < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2) {
		Val[k] = Val[0];
		Val[k + 1] = 0.;
		}
		}
		/*
		Message::Info("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++ ;		Num_Equ++;
} /* for i_Loop ... */		} /* for i_Loop ... */
} /* for i_MCPR ... */		} /* for i_MCPR ... */
List_Delete(DofGlobal_Equ_L) ;		List_Delete(DofGlobal_Equ_L);
List_Delete(DofGlobal_DofNode_L) ; List_Delete(DofGlobal_DofLoop_L) ;		List_Delete(DofGlobal_DofNode_L);
List_Delete(RegionIndex_L) ;		List_Delete(DofGlobal_DofLoop_L);
		List_Delete(RegionIndex_L);
}		}
/* ------------------------------------------------------------------------ */		/* ------------------------------------------------------------------------ */
/* T r e a t m e n t _ F e m F o r m u l a t i o n */		/* 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)		void Treatment_FemFormulation(struct Formulation *Formulation_P)
{		{
struct Element Element ;		struct Element Element;
struct QuantityStorage * QuantityStorage_P0, * QuantityStorage_P ;		struct QuantityStorage *QuantityStorage_P0, *QuantityStorage_P;
struct QuantityStorage QuantityStorage_S, QuantityStorage_GlobalEqu_S ;		struct QuantityStorage QuantityStorage_S, QuantityStorage_GlobalEqu_S;
struct Dof DofForNoDof_P [NBR_MAX_HARMONIC] ;		struct Dof DofForNoDof_P[NBR_MAX_HARMONIC];
struct EquationTerm * EquationTerm_P0 , * EquationTerm_P ;		struct EquationTerm *EquationTerm_P0, *EquationTerm_P;
struct GlobalQuantity * GlobalQuantity_P ;		struct GlobalQuantity *GlobalQuantity_P;
int Nbr_DefineQuantity ;		int Nbr_DefineQuantity;
struct DefineQuantity * DefineQuantity_P0 , * DefineQuantity_P ;		struct DefineQuantity *DefineQuantity_P0, *DefineQuantity_P;
List_T * FemLocalTermActive_L ;		List_T *FemLocalTermActive_L;
struct FemLocalTermActive FemLocalTermActive_S ;		struct FemLocalTermActive FemLocalTermActive_S;
List_T * QuantityStorage_L ;		List_T *QuantityStorage_L;
struct FemGlobalTermActive FemGlobalTermActive_S;		struct FemGlobalTermActive FemGlobalTermActive_S;
struct Group * GroupIn_P ;		struct Group *GroupIn_P;
int i, j, Nbr_Element, i_Element, Nbr_EquationTerm, i_EquTerm ;		int i, j, Nbr_Element, i_Element, Nbr_EquationTerm, i_EquTerm;
int Index_DefineQuantity, TraceGroupIndex_DefineQuantity ;		int Index_DefineQuantity, TraceGroupIndex_DefineQuantity;
List_T * InitialListInIndex_L ;		List_T *InitialListInIndex_L;
int Nbr_Region, i_Region, Num_Region ;		int Nbr_Region, i_Region, Num_Region;
int i_Region_Dof = 0, Num_Region_Dof = 0, i_Region_Dof_ini = 0;		int i_Region_Dof = 0, Num_Region_Dof = 0, i_Region_Dof_ini = 0;
int i_Region_Dof_end = 0, i_Region_Dof_skip = 0;		int i_Region_Dof_end = 0, i_Region_Dof_skip = 0;
extern struct Group * Generate_Group ;		extern struct Group *Generate_Group;
extern double ** MH_Moving_Matrix ;		extern double **MH_Moving_Matrix;
gMatrix A = Current.DofData->A ;		gMatrix A = Current.DofData->A;
gVector b = Current.DofData->b ;		gVector b = Current.DofData->b;
int Flag_Only ;		int Flag_Only;
/* --------------------------------------------------------------- */		/* --------------------------------------------------------------- */
/* 0. Initialization of an active zone (QuantityStorage) for each */		/* 0. Initialization of an active zone (QuantityStorage) for each */
/* DefineQuantity */		/* DefineQuantity */
/* --------------------------------------------------------------- */		/* --------------------------------------------------------------- */
if (!(Nbr_EquationTerm = List_Nbr(Formulation_P->Equation))){		if(!(Nbr_EquationTerm = List_Nbr(Formulation_P->Equation))) {
Message::Error("No equation in Formulation '%s'", Formulation_P->Name);		Message::Error("No equation in Formulation '%s'", Formulation_P->Name);
return;		return;
}		}
if (!(Nbr_DefineQuantity = List_Nbr(Formulation_P->DefineQuantity))){		if(!(Nbr_DefineQuantity = List_Nbr(Formulation_P->DefineQuantity))) {
Message::Error("No Quantity in Formulation '%s'", Formulation_P->Name);		Message::Error("No Quantity in Formulation '%s'", Formulation_P->Name);
return;		return;
}		}
DefineQuantity_P0 = (struct DefineQuantity*)		DefineQuantity_P0 =
List_Pointer(Formulation_P->DefineQuantity, 0) ;		(struct DefineQuantity *)List_Pointer(Formulation_P->DefineQuantity, 0);
QuantityStorage_L = List_Create(Nbr_DefineQuantity, 1,		QuantityStorage_L =
sizeof (struct QuantityStorage) ) ;		List_Create(Nbr_DefineQuantity, 1, sizeof(struct QuantityStorage));
QuantityStorage_S.NumLastElementForFunctionSpace =		QuantityStorage_S.NumLastElementForFunctionSpace =
QuantityStorage_S.NumLastElementForDofDefinition =		QuantityStorage_S.NumLastElementForDofDefinition =
QuantityStorage_S.NumLastElementForEquDefinition = 0 ;		QuantityStorage_S.NumLastElementForEquDefinition = 0;
for (i = 0 ; i < Nbr_DefineQuantity ; i++) {		for(i = 0; i < Nbr_DefineQuantity; i++) {
QuantityStorage_S.DefineQuantity = DefineQuantity_P0 + i ;		QuantityStorage_S.DefineQuantity = DefineQuantity_P0 + i;
if(QuantityStorage_S.DefineQuantity->Type == INTEGRALQUANTITY &&		if(QuantityStorage_S.DefineQuantity->Type == INTEGRALQUANTITY &&
QuantityStorage_S.DefineQuantity->IntegralQuantity.DefineQuantityIndexDof		QuantityStorage_S.DefineQuantity->IntegralQuantity
< 0){		.DefineQuantityIndexDof < 0) {
QuantityStorage_S.FunctionSpace = NULL ;		QuantityStorage_S.FunctionSpace = NULL;
QuantityStorage_S.TypeQuantity = VECTOR ; /* to change */		QuantityStorage_S.TypeQuantity = VECTOR; /* to change */
}		}
else{		else {
QuantityStorage_S.FunctionSpace = (struct FunctionSpace*)		QuantityStorage_S.FunctionSpace = (struct FunctionSpace *)List_Pointer(
List_Pointer(Problem_S.FunctionSpace,		Problem_S.FunctionSpace,
QuantityStorage_S.DefineQuantity->FunctionSpaceIndex) ;		QuantityStorage_S.DefineQuantity->FunctionSpaceIndex);
QuantityStorage_S.TypeQuantity = QuantityStorage_S.FunctionSpace->Type ;		QuantityStorage_S.TypeQuantity = QuantityStorage_S.FunctionSpace->Type;
}		}
List_Add(QuantityStorage_L, &QuantityStorage_S) ;		List_Add(QuantityStorage_L, &QuantityStorage_S);
}		}
QuantityStorage_P0 = (struct QuantityStorage*)List_Pointer(QuantityStorage_L,		QuantityStorage_P0 =
0) ;		(struct QuantityStorage *)List_Pointer(QuantityStorage_L, 0);
Get_InitDofOfElement(&Element) ;		Get_InitDofOfElement(&Element);
/* --------------------------------------------------------------- */		/* --------------------------------------------------------------- */
/* 1. Initialization of equation terms */		/* 1. Initialization of equation terms */
/* --------------------------------------------------------------- */		/* --------------------------------------------------------------- */
EquationTerm_P0 = (struct EquationTerm*)List_Pointer(Formulation_P->Equation,		EquationTerm_P0 =
0) ;		(struct EquationTerm *)List_Pointer(Formulation_P->Equation, 0);
FemLocalTermActive_L = List_Create(Nbr_EquationTerm, 1,		FemLocalTermActive_L =
sizeof (struct FemLocalTermActive) ) ;		List_Create(Nbr_EquationTerm, 1, sizeof(struct FemLocalTermActive));
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {		for(i_EquTerm = 0; i_EquTerm < Nbr_EquationTerm; i_EquTerm++) {
List_Add(FemLocalTermActive_L, &FemLocalTermActive_S) ;		List_Add(FemLocalTermActive_L, &FemLocalTermActive_S);
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;		EquationTerm_P = EquationTerm_P0 + i_EquTerm;
switch(EquationTerm_P->Type){		switch(EquationTerm_P->Type) {
		case GALERKIN:
case GALERKIN :		EquationTerm_P->Case.LocalTerm.Active =
EquationTerm_P->Case.LocalTerm.Active = (struct FemLocalTermActive*)		(struct FemLocalTermActive *)List_Pointer(FemLocalTermActive_L,
List_Pointer(FemLocalTermActive_L, i_EquTerm) ;		i_EquTerm);
switch (TreatmentStatus) {		switch(TreatmentStatus) {
case STATUS_PRE :		case STATUS_PRE:
Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;		Pre_InitTermOfFemEquation(EquationTerm_P, QuantityStorage_P0);
break ;		break;
case STATUS_CAL :		case STATUS_CAL:
Cal_InitGalerkinTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,		Cal_InitGalerkinTermOfFemEquation(EquationTerm_P, QuantityStorage_P0,
&QuantityStorage_S, DofForNoDof_P) ;		&QuantityStorage_S, DofForNoDof_P);
break ;		break;
}		}
break;		break;
case GLOBALTERM :		case GLOBALTERM:
switch (TreatmentStatus) {		switch(TreatmentStatus) {
case STATUS_PRE :		case STATUS_PRE:
Pre_InitGlobalTermOfFemEquation(EquationTerm_P, QuantityStorage_P0) ;		Pre_InitGlobalTermOfFemEquation(EquationTerm_P, QuantityStorage_P0);
break ;		break;
}		}
break;		break;
case GLOBALEQUATION :		case GLOBALEQUATION: break;
break ;
default :		default: Message::Error("Unknown type of equation term"); break;
Message::Error("Unknown type of equation term") ;
break ;
}		}
}		}
/* ---------------------------------------------------------- */		/* ---------------------------------------------------------- */
/* 2. Loop on geometrical elements : */		/* 2. Loop on geometrical elements : */
/* Treatment of eventual GALERKIN terms */		/* Treatment of eventual GALERKIN terms */
/* --------------------------------------------------------- */		/* --------------------------------------------------------- */
Nbr_Element = Geo_GetNbrGeoElements() ;		Nbr_Element = Geo_GetNbrGeoElements();
Message::ResetProgressMeter();		Message::ResetProgressMeter();
for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {		for(i_Element = 0; i_Element < Nbr_Element; i_Element++) {
		if(Generate_Group) {
if (Generate_Group) {		Element.Region = Geo_GetGeoElement(i_Element)->Region;
Element.Region = Geo_GetGeoElement(i_Element)->Region ;		while(
while (i_Element < Nbr_Element &&		i_Element < Nbr_Element &&
!List_Search(Generate_Group->InitialList,		!List_Search(Generate_Group->InitialList, &Element.Region, fcmp_int)) {
&Element.Region, fcmp_int) ) {		i_Element++;
i_Element++ ;		if(i_Element < Nbr_Element)
if (i_Element < Nbr_Element) Element.Region = Geo_GetGeoElement(i_Element		Element.Region = Geo_GetGeoElement(i_Element)->Region;
)->Region ;
}		}
if (i_Element == Nbr_Element) break ;		if(i_Element == Nbr_Element) break;
}		}
Element.GeoElement = Geo_GetGeoElement(i_Element) ;		Element.GeoElement = Geo_GetGeoElement(i_Element);
Element.Num = Element.GeoElement->Num ;		Element.Num = Element.GeoElement->Num;
Element.Type = Element.GeoElement->Type ;		Element.Type = Element.GeoElement->Type;
Current.Region = Element.Region = Element.GeoElement->Region ;		Current.Region = Element.Region = Element.GeoElement->Region;
/* ---------------------------- */		/* ---------------------------- */
/* 2.1. Loop on equation terms */		/* 2.1. Loop on equation terms */
/* ---------------------------- */		/* ---------------------------- */
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {		for(i_EquTerm = 0; i_EquTerm < Nbr_EquationTerm; i_EquTerm++) {
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;		EquationTerm_P = EquationTerm_P0 + i_EquTerm;
if (EquationTerm_P->Type == GALERKIN) {
/* if the element is in the support of integration of the term */		if(EquationTerm_P->Type == GALERKIN) {
/* ----------------------------------------------------------- */		/* if the element is in the support of integration of the term */
		/* ----------------------------------------------------------- */
GroupIn_P = (struct Group *)
List_Pointer(Problem_S.Group,		GroupIn_P = (struct Group *)List_Pointer(
EquationTerm_P->Case.LocalTerm.InIndex);		Problem_S.Group, EquationTerm_P->Case.LocalTerm.InIndex);
if ((GroupIn_P->Type != ELEMENTLIST &&		if((GroupIn_P->Type != ELEMENTLIST &&
List_Search(GroupIn_P->InitialList, &Element.Region, fcmp_int))		List_Search(GroupIn_P->InitialList, &Element.Region, fcmp_int)) ||
||		(GroupIn_P->Type == ELEMENTLIST &&
(GroupIn_P->Type == ELEMENTLIST &&		Check_IsEntityInExtendedGroup(GroupIn_P, Element.Num, 0))) {
Check_IsEntityInExtendedGroup(GroupIn_P, Element.Num, 0))		if(Message::GetVerbosity() == 10)
) {		printf("==> Element #%d, EquationTerm #%d/%d\n", Element.Num,
		i_EquTerm + 1, Nbr_EquationTerm);
/*
if (List_Search(((struct Group *)		Current.IntegrationSupportIndex =
List_Pointer(Problem_S.Group,		EquationTerm_P->Case.LocalTerm.InIndex;
EquationTerm_P->Case.
LocalTerm.InIndex))->InitialList,		if(EquationTerm_P->Case.LocalTerm.SubRegion >= 0) {
&Element.Region, fcmp_int ) ) {		struct Group *GroupSubRegion_P = (struct Group *)List_Pointer(
*/		Problem_S.Group, EquationTerm_P->Case.LocalTerm.SubRegion);
		if(List_Nbr(GroupSubRegion_P->InitialList) == 1) {
if(Message::GetVerbosity() == 10)		List_Read(GroupSubRegion_P->InitialList, 0, &Current.SubRegion);
printf("==> Element #%d, EquationTerm #%d/%d\n",
Element.Num, i_EquTerm+1, Nbr_EquationTerm) ;
Current.IntegrationSupportIndex = EquationTerm_P->Case.LocalTerm.InInde
x ;
if (EquationTerm_P->Case.LocalTerm.SubRegion >=0) {
struct Group * GroupSubRegion_P = (struct Group *)
List_Pointer(Problem_S.Group,
EquationTerm_P->Case.LocalTerm.SubRegion);
if (List_Nbr(GroupSubRegion_P->InitialList) == 1) {
List_Read(GroupSubRegion_P->InitialList, 0, &Current.SubRegion) ;
}		}
else {		else {
Message::Error("One region allowed in SubRegion");		Message::Error("One region allowed in SubRegion");
Current.SubRegion = -1;		Current.SubRegion = -1;
}		}
}		}
else		else
Current.SubRegion = -1;		Current.SubRegion = -1;
/* ---------------------------------------------------------- */		/* ---------------------------------------------------------- */
/* 2.1.1. Loop on quantities (test fcts and shape functions) */		/* 2.1.1. Loop on quantities (test fcts and shape functions) */
/* ---------------------------------------------------------- */		/* ---------------------------------------------------------- */
for (i = 0 ; i < EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex ;		std::vector<struct QuantityStorage *> QuantityStorageTrace;
i++) {
		for(i = 0; i < EquationTerm_P->Case.LocalTerm.Term.NbrQuantityIndex;
Index_DefineQuantity =		i++) {
EquationTerm_P->Case.LocalTerm.Term.QuantityIndexTable[i] ;		Index_DefineQuantity =
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;		EquationTerm_P->Case.LocalTerm.Term.QuantityIndexTable[i];
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;		DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity;
		QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity;
TraceGroupIndex_DefineQuantity =
EquationTerm_P->Case.LocalTerm.Term.QuantityTraceGroupIndexTable[i]		TraceGroupIndex_DefineQuantity = EquationTerm_P->Case.LocalTerm.Term
;		.QuantityTraceGroupIndexTable[i];
/* Only one analysis for each function space */		/* Only one analysis for each function space, unless we have a Trace
		(note that a quantity involved in a Trace cannot appear in
/*		the same term outside of a Trace) */
* Attention : l'operateur de trace ne fonctionne que si le champ
* dont on prend la trace n'intervient qu'une seule fois dans le term		if(QuantityStorage_P->NumLastElementForFunctionSpace !=
e.		Element.Num ||
* du a - manque de generalite du code au niveau de la gestion des		TraceGroupIndex_DefineQuantity >= 0) {
* espaces fonctionnels des fcts tests pour 'Trace de Dof'		switch(DefineQuantity_P->Type) {
* - et Christophe fatigue		case LOCALQUANTITY:
*/		if(TraceGroupIndex_DefineQuantity >= 0 &&
		Get_InitElementTrace(&Element,
if (QuantityStorage_P->NumLastElementForFunctionSpace != Element.Num		TraceGroupIndex_DefineQuantity)) {
||		/* Fill QuantityStorage for the first trace element; others
TraceGroupIndex_DefineQuantity >= 0 ) {		(e.g. for mortaring on nonconformal meshes) will be handled
		in a loop over the Cal_GalerkinTermOfFemEquation */
QuantityStorage_P->NumLastElementForFunctionSpace = Element.Num ;		QuantityStorage_P->NumLastElementForFunctionSpace =
		Element.ElementTrace->Num;
switch (DefineQuantity_P->Type) {		Get_DofOfElement(
case LOCALQUANTITY :		Element.ElementTrace, QuantityStorage_P->FunctionSpace,
if(TraceGroupIndex_DefineQuantity >= 0){		QuantityStorage_P, DefineQuantity_P->IndexInFunctionSpace);
Get_ElementTrace(&Element, TraceGroupIndex_DefineQuantity) ;		QuantityStorageTrace.push_back(QuantityStorage_P);
QuantityStorage_P->NumLastElementForFunctionSpace = Element.Ele		}
mentTrace->Num ;		else {
Get_DofOfElement		QuantityStorage_P->NumLastElementForFunctionSpace =
(Element.ElementTrace, QuantityStorage_P->FunctionSpace, Quan		Element.Num;
tityStorage_P,		Get_DofOfElement(&Element, QuantityStorage_P->FunctionSpace,
DefineQuantity_P->IndexInFunctionSpace) ;		QuantityStorage_P,
}		DefineQuantity_P->IndexInFunctionSpace);
else{		}
Get_DofOfElement		break;
(&Element, QuantityStorage_P->FunctionSpace, QuantityStorage_		case INTEGRALQUANTITY:
P,		QuantityStorage_P->NbrElementaryBasisFunction = 0;
DefineQuantity_P->IndexInFunctionSpace) ;		break;
}		default:
break ;		Message::Error("Bad kind of Quantity in Formulation '%s'",
case INTEGRALQUANTITY :
QuantityStorage_P->NbrElementaryBasisFunction = 0 ;
break ;
default :
Message::Error("Bad kind of Quantity in Formulation '%s'",
Formulation_P->Name);		Formulation_P->Name);
break;		break;
}		}
}		}
} /* for i = 0, 1 ... */		} /* for i = 0, 1 ... */
/* -------------------------------------- */		/* -------------------------------------- */
/* 2.1.2. Treatment of the equation term */		/* 2.1.2. Treatment of the equation term */
/* -------------------------------------- */		/* -------------------------------------- */
switch (TreatmentStatus) {		switch(TreatmentStatus) {
case STATUS_PRE :		case STATUS_PRE:
Pre_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0) ;		Pre_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0);
break ;		break;
case STATUS_CAL :		case STATUS_CAL:
Flag_Only = 0 ;		Flag_Only = 0;
if (Current.DofData->Flag_Only){		if(Current.DofData->Flag_Only) {
A = Current.DofData->A ;		A = Current.DofData->A;
b = Current.DofData->b ;		b = Current.DofData->b;
if (EquationTerm_P->Case.LocalTerm.MatrixIndex == -1)		if(EquationTerm_P->Case.LocalTerm.MatrixIndex == -1)
EquationTerm_P->Case.LocalTerm.MatrixIndex = 0 ;		EquationTerm_P->Case.LocalTerm.MatrixIndex = 0;
j = List_ISearch(Current.DofData->OnlyTheseMatrices,		j = List_ISearch(Current.DofData->OnlyTheseMatrices,
&EquationTerm_P->Case.LocalTerm.MatrixIndex, fcmp_		&EquationTerm_P->Case.LocalTerm.MatrixIndex,
int) ;		fcmp_int);
if(j!=-1){		if(j != -1) {
Flag_Only = 1 ;		Flag_Only = 1;
switch(EquationTerm_P->Case.LocalTerm.MatrixIndex){		switch(EquationTerm_P->Case.LocalTerm.MatrixIndex) {
case 1 :		case 1:
Current.DofData->A = Current.DofData->A1 ;		Current.DofData->A = Current.DofData->A1;
Current.DofData->b = Current.DofData->b1 ;		Current.DofData->b = Current.DofData->b1;
break;		break;
case 2 :		case 2:
Current.DofData->A = Current.DofData->A2 ;		Current.DofData->A = Current.DofData->A2;
Current.DofData->b = Current.DofData->b2 ;		Current.DofData->b = Current.DofData->b2;
break;		break;
case 3 :		case 3:
Current.DofData->A = Current.DofData->A3 ;		Current.DofData->A = Current.DofData->A3;
Current.DofData->b = Current.DofData->b3 ;		Current.DofData->b = Current.DofData->b3;
break;		break;
}		}
}		}
}/* Only the matrices that vary are recalculated */		} /* Only the matrices that vary are recalculated */
if (!Current.DofData->Flag_Only || (Current.DofData->Flag_Only && Fla
g_Only) ){
QuantityStorage_P = QuantityStorage_P0 +
EquationTerm_P->Case.LocalTerm.Term.DefineQuantityIndexEqu ;
if(EquationTerm_P->Type == GALERKIN){		if(!Current.DofData->Flag_Only ||
		(Current.DofData->Flag_Only && Flag_Only)) {
		if(EquationTerm_P->Type == GALERKIN) {
#if defined(HAVE_SMALLFEM)		#if defined(HAVE_SMALLFEM)
Cal_SmallFemTermOfFemEquation(&Element, EquationTerm_P, Quantity		Cal_SmallFemTermOfFemEquation(&Element, EquationTerm_P,
Storage_P0) ;		QuantityStorage_P0);
#else		#else
Cal_GalerkinTermOfFemEquation(&Element, EquationTerm_P, Quantity		Cal_GalerkinTermOfFemEquation(&Element, EquationTerm_P,
Storage_P0) ;		QuantityStorage_P0);
		/* if multiple candidate Trace elements, assemble the additional
		* contributions */
		while(Get_NextElementTrace(&Element)) {
		for(std::size_t k = 0; k < QuantityStorageTrace.size(); k++) {
		QuantityStorage_P = QuantityStorageTrace[k];
		QuantityStorage_P->NumLastElementForFunctionSpace =
		Element.ElementTrace->Num;
		Get_DofOfElement(Element.ElementTrace,
		QuantityStorage_P->FunctionSpace,
		QuantityStorage_P,
		DefineQuantity_P->IndexInFunctionSpace);
		}
		Cal_GalerkinTermOfFemEquation(&Element, EquationTerm_P,
		QuantityStorage_P0);
		}
#endif		#endif
}		}
if (Current.DofData->Flag_Only && Flag_Only){		if(Current.DofData->Flag_Only && Flag_Only) {
Current.DofData->A = A ;		Current.DofData->A = A;
Current.DofData->b = b ;		Current.DofData->b = b;
}		}
}/* Flag_Only */
break ;
case STATUS_CST :		} /* Flag_Only */
Cst_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0)		break;
;
break ;
}		case STATUS_CST:
		Cst_TermOfFemEquation(&Element, EquationTerm_P, QuantityStorage_P0);
		break;
		}
}/* if Support ... */		} /* if Support ... */
} /* if GALERKIN ... */		} /* if GALERKIN ... */
} /* for i_EquTerm ... */		} /* for i_EquTerm ... */
Message::ProgressMeter(i_Element + 1, Nbr_Element, (TreatmentStatus == STATU		Message::ProgressMeter(i_Element + 1, Nbr_Element,
S_PRE) ?		(TreatmentStatus == STATUS_PRE) ?
"Pre-processing" : "Processing (Generate)");		"Pre-processing" :
		"Processing (Generate)");
if(Message::GetErrorCount()) break;		if(Message::GetErrorCount()) break;
} /* for i_Element ... */		} /* for i_Element ... */
if (MH_Moving_Matrix) {		if(MH_Moving_Matrix) {
List_Delete(FemLocalTermActive_L) ; List_Delete(QuantityStorage_L) ;		List_Delete(FemLocalTermActive_L);
		List_Delete(QuantityStorage_L);
return;		return;
}		}
/* ------------------------------------------------------ */		/* ------------------------------------------------------ */
/* 3. Loop on equation terms : */		/* 3. Loop on equation terms : */
/* Treatment of eventual GLOBAL terms */		/* Treatment of eventual GLOBAL terms */
/* ------------------------------------------------------ */		/* ------------------------------------------------------ */
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {		for(i_EquTerm = 0; i_EquTerm < Nbr_EquationTerm; i_EquTerm++) {
		EquationTerm_P = EquationTerm_P0 + i_EquTerm;
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;
if (EquationTerm_P->Type == GLOBALTERM) {
		if(EquationTerm_P->Type == GLOBALTERM) {
EquationTerm_P->Case.GlobalTerm.Active = &FemGlobalTermActive_S;		EquationTerm_P->Case.GlobalTerm.Active = &FemGlobalTermActive_S;
InitialListInIndex_L =		InitialListInIndex_L =
((struct Group *)List_Pointer(Problem_S.Group,		((struct Group *)List_Pointer(Problem_S.Group,
EquationTerm_P->Case.GlobalTerm.InIndex))		EquationTerm_P->Case.GlobalTerm.InIndex))
->InitialList ;		->InitialList;
List_Sort(InitialListInIndex_L, fcmp_int) ;		List_Sort(InitialListInIndex_L, fcmp_int);
Nbr_Region = List_Nbr(InitialListInIndex_L) ;		Nbr_Region = List_Nbr(InitialListInIndex_L);
/* ---------------------------------------------- */		/* ---------------------------------------------- */
/* 3.1. Loop on Regions belonging to the support */		/* 3.1. Loop on Regions belonging to the support */
/* ---------------------------------------------- */		/* ---------------------------------------------- */
for (i_Region = 0 ; i_Region < Nbr_Region ; i_Region++) {		for(i_Region = 0; i_Region < Nbr_Region; i_Region++) {
List_Read(InitialListInIndex_L, i_Region, &Num_Region) ;		List_Read(InitialListInIndex_L, i_Region, &Num_Region);
Current.Region = Num_Region ;		Current.Region = Num_Region;
switch (EquationTerm_P->Case.GlobalTerm.SubType) {		switch(EquationTerm_P->Case.GlobalTerm.SubType) {
case EQ_ST_SELF:		case EQ_ST_SELF:
i_Region_Dof_ini = i_Region; i_Region_Dof_end = i_Region+1;		i_Region_Dof_ini = i_Region;
		i_Region_Dof_end = i_Region + 1;
i_Region_Dof_skip = -1;		i_Region_Dof_skip = -1;
break;		break;
case EQ_ST_MUTUAL:		case EQ_ST_MUTUAL:
i_Region_Dof_ini = 0; i_Region_Dof_end = Nbr_Region;		i_Region_Dof_ini = 0;
		i_Region_Dof_end = Nbr_Region;
i_Region_Dof_skip = i_Region;		i_Region_Dof_skip = i_Region;
break;		break;
case EQ_ST_SELFANDMUTUAL:		case EQ_ST_SELFANDMUTUAL:
i_Region_Dof_ini = 0; i_Region_Dof_end = Nbr_Region;		i_Region_Dof_ini = 0;
		i_Region_Dof_end = Nbr_Region;
i_Region_Dof_skip = -1;		i_Region_Dof_skip = -1;
break;		break;
}		}
// Possible mutual terms (need of double-piece-wise Function fct[r1,r2])		// Possible mutual terms (need of double-piece-wise Function fct[r1,r2])
for (i_Region_Dof = i_Region_Dof_ini; i_Region_Dof < i_Region_Dof_end;		for(i_Region_Dof = i_Region_Dof_ini; i_Region_Dof < i_Region_Dof_end;
i_Region_Dof++) {		i_Region_Dof++) {
		if(i_Region_Dof != i_Region_Dof_skip) {
if (i_Region_Dof != i_Region_Dof_skip) {		List_Read(InitialListInIndex_L, i_Region_Dof, &Num_Region_Dof);
		Current.SubRegion =
List_Read(InitialListInIndex_L, i_Region_Dof, &Num_Region_Dof) ;		Num_Region_Dof; // used in double-piece-wise Function
Current.SubRegion = Num_Region_Dof; // used in double-piece-wise Fun
ction		/* ----------------------------------------------------------------
		*/
/* ----------------------------------------------------------------		/* 3.1.1. Loop on Quantities (test functions and shape functions)
*/		*/
/* 3.1.1. Loop on Quantities (test functions and shape functions)		/* ----------------------------------------------------------------
*/		*/
/* ----------------------------------------------------------------
*/		for(i = 0;
		i < EquationTerm_P->Case.GlobalTerm.Term.NbrQuantityIndex;
for (i = 0 ; i < EquationTerm_P->Case.GlobalTerm.Term.NbrQuantityInd		i++) {
ex ; i++) {
Index_DefineQuantity =		Index_DefineQuantity =
EquationTerm_P->Case.GlobalTerm.Term.QuantityIndexTable[i] ;		EquationTerm_P->Case.GlobalTerm.Term.QuantityIndexTable[i];
DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity ;		DefineQuantity_P = DefineQuantity_P0 + Index_DefineQuantity;
QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity ;		QuantityStorage_P = QuantityStorage_P0 + Index_DefineQuantity;
GlobalQuantity_P = (struct GlobalQuantity*)		GlobalQuantity_P = (struct GlobalQuantity *)List_Pointer(
List_Pointer(QuantityStorage_P->FunctionSpace->GlobalQuantity,		QuantityStorage_P->FunctionSpace->GlobalQuantity,
*(int*)List_Pointer(DefineQuantity_P->IndexInFuncti		*(int *)List_Pointer(DefineQuantity_P->IndexInFunctionSpace,
onSpace, 0)) ;		0));
/* Only one Function space analysis */		/* Only one Function space analysis */
/* -------------------------------- */		/* -------------------------------- */
if (QuantityStorage_P->NumLastElementForFunctionSpace != Num_Regio		if(QuantityStorage_P->NumLastElementForFunctionSpace !=
n_Dof) {		Num_Region_Dof) {
QuantityStorage_P->NumLastElementForFunctionSpace = Num_Region_D		QuantityStorage_P->NumLastElementForFunctionSpace =
of ;		Num_Region_Dof;
switch (DefineQuantity_P->Type) {		switch(DefineQuantity_P->Type) {
case GLOBALQUANTITY :		case GLOBALQUANTITY:
Get_DofOfRegion		Get_DofOfRegion(Num_Region_Dof, GlobalQuantity_P,
(Num_Region_Dof, GlobalQuantity_P,		QuantityStorage_P->FunctionSpace,
QuantityStorage_P->FunctionSpace, QuantityStorage_P) ;		QuantityStorage_P);
break ;		break;
default :		default:
Message::Error("Bad kind of Quantity in Formulation '%s'",		Message::Error("Bad kind of Quantity in Formulation '%s'",
Formulation_P->Name);		Formulation_P->Name);
break;		break;
}		}
}		}
// Particular QuantityStorage for Equ		// Particular QuantityStorage for Equ
if (Num_Region != Num_Region_Dof &&		if(Num_Region != Num_Region_Dof &&
Index_DefineQuantity ==		Index_DefineQuantity == EquationTerm_P->Case.GlobalTerm.Term
EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexEqu) {		.DefineQuantityIndexEqu) {
switch (DefineQuantity_P->Type) {		switch(DefineQuantity_P->Type) {
case GLOBALQUANTITY :		case GLOBALQUANTITY:
Get_DofOfRegion		Get_DofOfRegion(Num_Region, GlobalQuantity_P,
(Num_Region, GlobalQuantity_P,		QuantityStorage_P->FunctionSpace,
QuantityStorage_P->FunctionSpace, &QuantityStorage_GlobalEq		&QuantityStorage_GlobalEqu_S);
u_S) ;		break;
break ;		default:
default :
Message::Error("Bad kind of Quantity in Formulation '%s'",		Message::Error("Bad kind of Quantity in Formulation '%s'",
Formulation_P->Name);		Formulation_P->Name);
break;		break;
}		}
}		}
} /* for i = 0, 1 ... */		} /* for i = 0, 1 ... */
// QuantityStorage for Equ and Dof (can differ for SubType Mutual)		// QuantityStorage for Equ and Dof (can differ for SubType Mutual)
EquationTerm_P->Case.GlobalTerm.Active->QuantityStorageEqu_P =		EquationTerm_P->Case.GlobalTerm.Active->QuantityStorageEqu_P =
(Num_Region == Num_Region_Dof)?		(Num_Region == Num_Region_Dof) ?
QuantityStorage_P0 +		QuantityStorage_P0 +
EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexEqu		EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexEqu :
:		&QuantityStorage_GlobalEqu_S;
&QuantityStorage_GlobalEqu_S
;
EquationTerm_P->Case.GlobalTerm.Active->flag_Dof =		EquationTerm_P->Case.GlobalTerm.Active->flag_Dof =
(EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexDof >= 0)		(EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexDof >=
;		0);
EquationTerm_P->Case.GlobalTerm.Active->QuantityStorageDof_P =		EquationTerm_P->Case.GlobalTerm.Active->QuantityStorageDof_P =
(EquationTerm_P->Case.GlobalTerm.Active->flag_Dof)?		(EquationTerm_P->Case.GlobalTerm.Active->flag_Dof) ?
QuantityStorage_P0 +		QuantityStorage_P0 +
EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexDof		EquationTerm_P->Case.GlobalTerm.Term.DefineQuantityIndexDof :
: NULL;		NULL;
/* ------------------------------ */		/* ------------------------------ */
/* 3.1.2. Treatment of the term */		/* 3.1.2. Treatment of the term */
/* ------------------------------ */		/* ------------------------------ */
switch (TreatmentStatus) {		switch(TreatmentStatus) {
case STATUS_PRE :		case STATUS_PRE:
Pre_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,		Pre_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
QuantityStorage_P0) ;		QuantityStorage_P0);
break ;		break;
case STATUS_CAL :		case STATUS_CAL:
Cal_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,		Cal_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
QuantityStorage_P0,		QuantityStorage_P0,
&QuantityStorage_S, DofForNoDof_P) ;		&QuantityStorage_S, DofForNoDof_P);
break ;		break;
case STATUS_CST :		case STATUS_CST:
Cst_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,		Cst_GlobalTermOfFemEquation(Num_Region, EquationTerm_P,
QuantityStorage_P0) ;		QuantityStorage_P0);
break ;		break;
}		}
}		}
} // for i_Region_Dof		} // for i_Region_Dof
} // for i_Region		} // for i_Region
} /* if GLOBALTERM ... */		} /* if GLOBALTERM ... */
} /* for i_EquTerm ... */		} /* for i_EquTerm ... */
/* --------------------------------------------------------- */		/* --------------------------------------------------------- */
/* 4. Loop on equation terms : */		/* 4. Loop on equation terms : */
/* Treatment of eventual GLOBAL EQUATION terms */		/* Treatment of eventual GLOBAL EQUATION terms */
/* --------------------------------------------------------- */		/* --------------------------------------------------------- */
for (i_EquTerm = 0 ; i_EquTerm < Nbr_EquationTerm ; i_EquTerm++) {		for(i_EquTerm = 0; i_EquTerm < Nbr_EquationTerm; i_EquTerm++) {
EquationTerm_P = EquationTerm_P0 + i_EquTerm ;		EquationTerm_P = EquationTerm_P0 + i_EquTerm;
if (EquationTerm_P->Type == GLOBALEQUATION) {		if(EquationTerm_P->Type == GLOBALEQUATION) {
		if(EquationTerm_P->Case.GlobalEquation.ConstraintIndex >= 0)
if (EquationTerm_P->Case.GlobalEquation.ConstraintIndex >= 0)		switch(TreatmentStatus) {
		case STATUS_PRE:
switch (TreatmentStatus) {		Pre_FemGlobalEquation(EquationTerm_P, DefineQuantity_P0,
case STATUS_PRE :		QuantityStorage_P0);
Pre_FemGlobalEquation(EquationTerm_P,		break;
DefineQuantity_P0, QuantityStorage_P0) ;		case STATUS_CAL:
break ;		Cal_FemGlobalEquation(EquationTerm_P, DefineQuantity_P0,
case STATUS_CAL :		QuantityStorage_P0);
Cal_FemGlobalEquation(EquationTerm_P,		break;
DefineQuantity_P0, QuantityStorage_P0) ;		}
break ;
}
} /* if GLOBALEQUATION ... */		} /* if GLOBALEQUATION ... */
} /* for i_EquTerm ... */		} /* for i_EquTerm ... */
/* -------------------------- */		/* -------------------------- */
/* 5. End of FEM treatment */		/* 5. End of FEM treatment */
/* -------------------------- */		/* -------------------------- */
List_Delete(FemLocalTermActive_L) ; List_Delete(QuantityStorage_L) ;		List_Delete(FemLocalTermActive_L);
		List_Delete(QuantityStorage_L);
Cal_EndGalerkinTermOfFemEquation();		Cal_EndGalerkinTermOfFemEquation();
}		}
/* ------------------------------------------------------------------------ */		/* ------------------------------------------------------------------------ */
/* 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 */		/* 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)		void Treatment_GlobalFormulation(struct Formulation *Formulation_P)
{		{
Message::Error("You should not be here!") ;		Message::Error("You should not be here!");
}		}
/* ------------------------------------------------------------------------ */		/* ------------------------------------------------------------------------ */
/* T r e a t m e n t _ F o r m u l a t i o n */		/* 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)		void Treatment_Formulation(struct Formulation *Formulation_P)
{		{
switch (Formulation_P->Type) {		switch(Formulation_P->Type) {
		case FEMEQUATION: Treatment_FemFormulation(Formulation_P); break;
		case GLOBALEQUATION: Treatment_GlobalFormulation(Formulation_P); break;
case FEMEQUATION :		default:
Treatment_FemFormulation(Formulation_P) ;		Message::Error("Unknown type for Formulation '%s'", Formulation_P->Name);
break ;		break;
case GLOBALEQUATION :
Treatment_GlobalFormulation(Formulation_P) ;
break ;
default :
Message::Error("Unknown type for Formulation '%s'", Formulation_P->Name) ;
break ;
}		}
}		}
End of changes. 106 change blocks.
601 lines changed or deleted		562 lines changed or added

---