#define RCSID "$Id: F_MultiHar.c,v 1.26 2006/02/26 00:42:53 geuzaine Exp $"
/*
* Copyright (C) 1997-2006 P. Dular, C. Geuzaine
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
* USA.
*
* Please report all bugs and problems to <getdp@geuz.org>.
*
* Contributor(s):
* Johan Gyselinck
*/
#include "GetDP.h"
#include "Data_DefineE.h"
#include "F_Function.h"
#include "GeoData.h"
#include "Get_Geometry.h"
#include "Cal_Value.h"
#include "Treatment_Formulation.h"
#include "Cal_Quantity.h"
#include "CurrentData.h"
#include "Numeric.h"
#include "Tools.h"
struct MH_InitData{
int Case ;
int NbrPoints, NbrPointsX; /* number of samples per smallest and fundamental period resp. */
struct DofData *DofData ;
double **H, ***HH ;
double *t, *w;
};
List_T * MH_InitData_L = NULL;
int fcmp_MH_InitData(const void * a, const void * b) {
int Result ;
if ((Result = ((struct MH_InitData *)a)->DofData - ((struct MH_InitData *)b)->DofData) != 0)
return Result ;
if ((Result = ((struct MH_InitData *)a)->Case - ((struct MH_InitData *)b)->Case) != 0)
return Result ;
if (((struct MH_InitData *)a)->Case != 3)
return ((struct MH_InitData *)a)->NbrPoints - ((struct MH_InitData *)b)->NbrPoints ;
else
return ((struct MH_InitData *)a)->NbrPointsX - ((struct MH_InitData *)b)->NbrPointsX ;
}
int NbrValues_Type (int Type){
switch (Type){
case SCALAR : return 1 ;
case VECTOR : case TENSOR_DIAG : return 3 ;
case TENSOR_SYM : return 6 ;
case TENSOR : return 9 ;
default :
Msg(GERROR, "Unknown type in NbrValues_Type");
return 0;
}
}
double Product_SCALARxSCALARxSCALAR (double *V1, double *V2, double *V3){
return V1[0] * V2[0] * V3[0] ;
}
double Product_VECTORxTENSOR_SYMxVECTOR (double *V1, double *V2, double *V3){
return V3[0] * (V1[0] * V2[0] + V1[1] * V2[1] + V1[2] * V2[2]) +
V3[1] * (V1[0] * V2[1] + V1[1] * V2[3] + V1[2] * V2[4]) +
V3[2] * (V1[0] * V2[2] + V1[1] * V2[4] + V1[2] * V2[5]) ;
}
double Product_VECTORxTENSOR_DIAGxVECTOR (double *V1, double *V2, double *V3){
return V1[0] * V2[0] * V3[0] + V1[1] * V2[1] * V3[1] + V1[2] * V2[2] * V3[2] ;
}
double Product_VECTORxSCALARxVECTOR (double *V1, double *V2, double *V3){
return V2[0] * (V1[0] * V3[0] + V1[1] * V3[1] + V1[2] * V3[2]) ;
}
void * Get_RealProductFunction_Type1xType2xType1 (int Type1, int Type2) {
/* Type1 * Type2 * Type1 */
GetDP_Begin("Get_RealProductFunction_Type1xType2xType1");
if (Type1 == SCALAR && Type2 == SCALAR) {
GetDP_Return((void *)Product_SCALARxSCALARxSCALAR) ;
} else if (Type1 == VECTOR && Type2 == TENSOR_SYM) {
GetDP_Return((void *)Product_VECTORxTENSOR_SYMxVECTOR) ;
} else if (Type1 == VECTOR && Type2 == TENSOR_DIAG) {
GetDP_Return((void *)Product_VECTORxTENSOR_DIAGxVECTOR) ;
} else if (Type1 == VECTOR && Type1 == SCALAR) {
GetDP_Return((void *)Product_VECTORxSCALARxVECTOR) ;
} else {
Msg(GERROR, "Not allowed types in Get_RealProductFunction_Type1xType2xType1");
GetDP_Return(NULL) ;
}
}
/* ------------------------------------------------------------------------ */
/* MH_Get_InitData */
/* ------------------------------------------------------------------------ */
/*
Case = 1 : MHTransform NbrPoints (samples per smallest period) is given,
NbrPointsX (samples per fundamental period) is derived
Case = 2 : MHJacNL NbrPoints given, NbrPointsX derived
Case = 3 : HarmonicToTime NbrPointsX given, NbrPoints derived
*/
void MH_Get_InitData(int Case, int NbrPoints, int *NbrPointsX_P,
double ***H_P, double ****HH_P, double **t_P, double **w_P){
int NbrHar, iPul, iTime, iHar, jHar, NbrPointsX ;
double *Val_Pulsation, MaxPuls, MinPuls ;
double **H, ***HH = 0, *t, *w ;
struct MH_InitData MH_InitData_S, *MH_InitData_P ;
GetDP_Begin("MH_Get_InitData(");
MH_InitData_S.Case = Case;
MH_InitData_S.DofData = Current.DofData;
MH_InitData_S.NbrPoints = NbrPoints;
MH_InitData_S.NbrPointsX = NbrPointsX = *NbrPointsX_P;
if (MH_InitData_L == NULL)
MH_InitData_L = List_Create(1, 1, sizeof(struct MH_InitData)) ;
if ((MH_InitData_P = (struct MH_InitData *)
List_PQuery(MH_InitData_L, &MH_InitData_S, fcmp_MH_InitData))){
*H_P = MH_InitData_P->H;
*HH_P = MH_InitData_P->HH;
*t_P = MH_InitData_P->t;
*w_P = MH_InitData_P->w;
*NbrPointsX_P = MH_InitData_P->NbrPointsX;
GetDP_End;
}
NbrHar = Current.NbrHar;
Val_Pulsation = Current.DofData->Val_Pulsation;
MaxPuls = 0. ; MinPuls = 1.e99 ;
for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) {
if (Val_Pulsation[iPul] && Val_Pulsation[iPul] < MinPuls)
MinPuls = Val_Pulsation[iPul] ;
if (Val_Pulsation[iPul] && Val_Pulsation[iPul] > MaxPuls)
MaxPuls = Val_Pulsation[iPul] ;
}
if (Case != 3)
NbrPointsX = (int)((MaxPuls/MinPuls*(double)NbrPoints));
else
NbrPoints = (int)((MinPuls/MaxPuls*(double)NbrPointsX));
Msg(INFO, "MH_Get_InitData => NbrHar = %d NbrPoints = %d|%d Case = %d",
NbrHar, NbrPoints, NbrPointsX, Case);
t = (double *)Malloc(sizeof(double)*NbrPointsX) ;
if (Case != 3)
for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
t[iTime] = (double)iTime/(double)NbrPointsX/(MinPuls/TWO_PI) ;
else
for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
t[iTime] = (double)iTime/((double)NbrPointsX-1.)/(MinPuls/TWO_PI) ;
/*
w = (double *)Malloc(sizeof(double)*NbrPointsX) ;
for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
w[iTime] = 2. / (double)NbrPointsX ;
*/
w = (double *)Malloc(sizeof(double)*NbrHar) ;
for (iPul = 0 ; iPul < NbrHar/2 ; iPul++)
if (Val_Pulsation[iPul]){
w[2*iPul ] = 2. / (double)NbrPointsX ;
w[2*iPul+1] = 2. / (double)NbrPointsX ;
}else{
w[2*iPul ] = 1. / (double)NbrPointsX ;
w[2*iPul+1] = 1. / (double)NbrPointsX ;
}
H = (double **)Malloc(sizeof(double *)*NbrPointsX) ;
for (iTime = 0 ; iTime < NbrPointsX ; iTime++){
H[iTime] = (double *)Malloc(sizeof(double)*NbrHar) ;
for (iPul = 0 ; iPul < NbrHar/2 ; iPul++) {
/* if (Val_Pulsation [iPul]){ */
H[iTime][2*iPul ] = cos(Val_Pulsation[iPul] * t[iTime]) ;
H[iTime][2*iPul+1] = - sin(Val_Pulsation[iPul] * t[iTime]) ;
}
/*
}
else {
H[iTime][2*iPul ] = 0.5 ;
H[iTime][2*iPul + 1] = 0 ;
}
*/
}
/*
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = iHar ; jHar < NbrHar ; jHar++){
sum = 0.;
for (iTime = 0 ; iTime < NbrPointsX ; iTime++)
sum += w[iTime] * H[iTime][iHar] * H[iTime][jHar] ;
sum -= (iHar==jHar)? 1. : 0. ;
printf("iHar %d jHar %d sum %e\n", iHar, jHar, sum);
}
*/
if (Case == 2) {
HH = (double ***)Malloc(sizeof(double **)*NbrPointsX) ;
for (iTime = 0 ; iTime < NbrPointsX ; iTime++){
HH[iTime] = (double **)Malloc(sizeof(double *)*NbrHar) ;
for (iHar = 0 ; iHar < NbrHar ; iHar++){
HH[iTime][iHar] = (double *)Malloc(sizeof(double)*NbrHar) ;
for (jHar = 0 ; jHar < NbrHar ; jHar++){
if (Val_Pulsation [iHar/2] && Val_Pulsation [jHar/2] )
HH[iTime][iHar][jHar] = 2. / (double)NbrPointsX * H[iTime][iHar] * H[iTime][jHar] ;
else
HH[iTime][iHar][jHar] = 1. / (double)NbrPointsX * H[iTime][iHar] * H[iTime][jHar] ;
}
}
}
}
*H_P = MH_InitData_S.H = H;
*t_P = MH_InitData_S.t = t;
*w_P = MH_InitData_S.w = w;
*HH_P = MH_InitData_S.HH = HH;
*NbrPointsX_P = MH_InitData_S.NbrPointsX = NbrPointsX;
List_Add (MH_InitData_L, &MH_InitData_S);
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* F_MHToTime0 (HarmonicToTime in PostOperation) */
/* ------------------------------------------------------------------------ */
void F_MHToTime0 (int init, struct Value * A, struct Value * V,
int iTime, int NbrPointsX, double * TimeMH) {
static double **H, ***HH, *t, *weight;
int iVal, nVal, iHar;
GetDP_Begin("F_MHToTime0");
if (Current.NbrHar == 1) GetDP_End ;
if (!init) MH_Get_InitData(3, 0, &NbrPointsX, &H, &HH, &t, &weight);
*TimeMH = t[iTime] ;
V->Type = A->Type ;
nVal = NbrValues_Type (A->Type) ;
for (iVal = 0 ; iVal < nVal ; iVal++){
V->Val[iVal] = 0;
for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
V->Val[iVal] += H[iTime][iHar] * A->Val[iHar*MAX_DIM+iVal] ;
}
}
/* ---------------------------------------------------------------------- */
/* F_MHToTime */
/* ---------------------------------------------------------------------- */
void F_MHToTime (struct Function * Fct, struct Value * A, struct Value * V) {
int iHar, iVal, nVal ;
double time, H[NBR_MAX_HARMONIC];
struct Value Vtemp;
GetDP_Begin("F_MHToTime");
if (Current.NbrHar == 1) Msg(GERROR, "'F_MHToTime' only for Multi-Harmonic stuff") ;
if((A+1)->Type != SCALAR)
Msg(GERROR, "'F_MHToTime' requires second scalar argument (time)");
time = (A+1)->Val[0] ;
for (iHar = 0 ; iHar < Current.NbrHar/2 ; iHar++) {
/* if (Current.DofData->Val_Pulsation [iHar]){ */
H[2*iHar ] = cos(Current.DofData->Val_Pulsation[iHar] * time) ;
H[2*iHar+1] = - sin(Current.DofData->Val_Pulsation[iHar] * time) ;
/*
}
else {
H[2*iHar ] = 0.5 ;
H[2*iHar+1] = 0 ;
}
*/
}
nVal = NbrValues_Type (A->Type) ;
for (iVal = 0 ; iVal < MAX_DIM ; iVal++)
for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
Vtemp.Val[iHar*MAX_DIM+iVal] = 0.;
for (iVal = 0 ; iVal < nVal ; iVal++)
for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
Vtemp.Val[iVal] += H[iHar] * A->Val[iHar*MAX_DIM+iVal] ;
V->Type = A->Type ;
for (iVal = 0 ; iVal < MAX_DIM ; iVal++)
for (iHar = 0 ; iHar < Current.NbrHar ; iHar++)
V->Val[iHar*MAX_DIM+iVal] = Vtemp.Val[iHar*MAX_DIM+iVal] ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* MHTransform */
/* ------------------------------------------------------------------------ */
void MHTransform(struct Element * Element, struct QuantityStorage * QuantityStorage_P0,
double u, double v, double w, struct Value *MH_Value,
struct Expression * Expression_P, int NbrPoints){
double **H, ***HH, *t, *weight ;
int NbrHar;
struct Value t_Value, MH_Value_Tr;
int NbrPointsX, iVal, nVal1, nVal2 = 0, iHar, iTime;
MH_Get_InitData(1, NbrPoints, &NbrPointsX, &H, &HH, &t, &weight);
nVal1 = NbrValues_Type (MH_Value->Type) ;
t_Value.Type = MH_Value_Tr.Type = MH_Value->Type;
NbrHar = Current.NbrHar; /* save NbrHar */
Current.NbrHar = 1; /* evaluation in time domain ! */
for (iVal = 0 ; iVal < MAX_DIM ; iVal++) for (iHar = 0 ; iHar < NbrHar ; iHar++)
MH_Value_Tr.Val[iHar*MAX_DIM+iVal] = 0. ;
for (iTime = 0 ; iTime < NbrPointsX ; iTime++) {
for (iVal = 0 ; iVal < nVal1 ; iVal++){ /* evaluation of MH_Value at iTime-th time point */
t_Value.Val[iVal] = 0;
for (iHar = 0 ; iHar < NbrHar ; iHar++)
t_Value.Val[iVal] += H[iTime][iHar] * MH_Value->Val[iHar*MAX_DIM+iVal] ;
}
/* evaluation of the function */
Get_ValueOfExpression(Expression_P, QuantityStorage_P0, u, v, w, &t_Value);
if (!iTime) nVal2 = NbrValues_Type (t_Value.Type) ;
for (iVal = 0 ; iVal < nVal2 ; iVal++) for (iHar = 0 ; iHar < NbrHar ; iHar++)
MH_Value_Tr.Val[iHar*MAX_DIM+iVal] +=
weight[iHar] * H[iTime][iHar] * t_Value.Val[iVal] ;
/* weight[iTime] * H[iTime][iHar] * t_Value.Val[iVal] ; */
} /* for iTime */
for (iVal = 0 ; iVal < nVal2 ; iVal++) for (iHar = 0 ; iHar < NbrHar ; iHar++)
MH_Value->Val[iHar*MAX_DIM+iVal] = MH_Value_Tr.Val[iHar*MAX_DIM+iVal] ;
MH_Value->Type = t_Value.Type ;
Current.NbrHar = NbrHar ;
GetDP_End ;
}
/* ----------------------------------------------------------------------------------- */
/* C a l _ I n i t G a l e r k i n T e r m O f F e m E q u a t i o n _ M H J a c N L */
/* ----------------------------------------------------------------------------------- */
void Cal_InitGalerkinTermOfFemEquation_MHJacNL(struct EquationTerm * EquationTerm_P) {
struct FemLocalTermActive * FI ;
List_T * WholeQuantity_L;
struct WholeQuantity *WholeQuantity_P0 ;
int i_WQ ;
GetDP_Begin("Cal_InitGalerkinTermOfFemEquation_MHJacNL");
FI = EquationTerm_P->Case.LocalTerm.Active ;
FI->MHJacNL = 0 ;
/* search for MHJacNL-term(s) */
WholeQuantity_L = EquationTerm_P->Case.LocalTerm.Term.WholeQuantity ;
WholeQuantity_P0 = (struct WholeQuantity*)List_Pointer(WholeQuantity_L, 0) ;
i_WQ = 0 ; while ( i_WQ < List_Nbr(WholeQuantity_L) &&
(WholeQuantity_P0 + i_WQ)->Type != WQ_MHJACNL) i_WQ++ ;
if (i_WQ == List_Nbr(WholeQuantity_L) )
GetDP_End ; /* no MHJacNL stuff, let's get the hell out of here ! */
/* check if Galerkin term produces symmetrical contribution to system matrix */
if (!FI->SymmetricalMatrix)
Msg(GERROR, "Galerkin term with MHJacNL must be symmetrical");
if (List_Nbr(WholeQuantity_L) == 3){
if (i_WQ != 0 ||
EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity != 1 ||
(WholeQuantity_P0 + 2)->Type != WQ_BINARYOPERATOR ||
(WholeQuantity_P0 + 2)->Case.Operator.TypeOperator != OP_TIME)
Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL");
FI->MHJacNL_Factor = 1.;
} else if (List_Nbr(WholeQuantity_L) == 5){
if ((WholeQuantity_P0 + 0)->Type != WQ_CONSTANT ||
i_WQ != 1 ||
(WholeQuantity_P0 + 2)->Type != WQ_BINARYOPERATOR ||
(WholeQuantity_P0 + 2)->Case.Operator.TypeOperator != OP_TIME ||
EquationTerm_P->Case.LocalTerm.Term.DofIndexInWholeQuantity != 3 ||
(WholeQuantity_P0 + 4)->Type != WQ_BINARYOPERATOR ||
(WholeQuantity_P0 + 4)->Case.Operator.TypeOperator != OP_TIME)
Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL");
FI->MHJacNL_Factor = WholeQuantity_P0->Case.Constant ;
/* printf(" Factor = %e \n" , FI->MHJacNL_Factor); */
} else {
Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL (%d terms) ",
List_Nbr(WholeQuantity_L));
}
if(EquationTerm_P->Case.LocalTerm.Term.CanonicalWholeQuantity_Equ != CWQ_NONE)
Msg(GERROR, "Not allowed expression in Galerkin term with MHJacNL");
if (EquationTerm_P->Case.LocalTerm.Term.TypeTimeDerivative != JACNL_)
Msg(GERROR, "MHJacNL can only be used with JACNL") ;
FI->MHJacNL = 1 ;
FI->MHJacNL_Index = (WholeQuantity_P0 + i_WQ)->Case.MHJacNL.Index ;
FI->MHJacNL_HarOffSet = 2 * (WholeQuantity_P0 + i_WQ)->Case.MHJacNL.FreqOffSet ;
if (FI->MHJacNL_HarOffSet > Current.NbrHar-2){
Msg(WARNING, "FreqOffSet in 'MHJacNL' cannot exceed %d => set to %d ",
Current.NbrHar/2-1, Current.NbrHar/2-1) ;
FI->MHJacNL_HarOffSet = Current.NbrHar-2 ;
}
MH_Get_InitData(2, (WholeQuantity_P0 + i_WQ)->Case.MHJacNL.NbrPoints,
&FI->MHJacNL_NbrPointsX, &FI->MHJacNL_H, &FI->MHJacNL_HH,
&FI->MHJacNL_t, &FI->MHJacNL_w) ;
GetDP_End ;
}
#define OFFSET (iHar < NbrHar-OffSet)? 0 : iHar-NbrHar+OffSet+2-iHar%2
/* --------------------------------------------------------------------------- */
/* C a l _ G a l e r k i n T e r m O f F e m E q u a t i o n _ M H J a c N L */
/* --------------------------------------------------------------------------- */
void Cal_GalerkinTermOfFemEquation_MHJacNL(struct Element * Element,
struct EquationTerm * EquationTerm_P,
struct QuantityStorage * QuantityStorage_P0) {
struct FemLocalTermActive * FI ;
struct QuantityStorage * QuantityStorage_P;
struct IntegrationCase * IntegrationCase_P ;
struct Quadrature * Quadrature_P ;
int Nbr_Dof, NbrHar ;
double vBFuDof [NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
double vBFxDof [NBR_MAX_BASISFUNCTIONS][MAX_DIM] ;
double Val_Dof [NBR_MAX_BASISFUNCTIONS][NBR_MAX_HARMONIC] ;
double Val [3*NBR_MAX_BASISFUNCTIONS];
int i, j, k, Type_Dimension, Nbr_IntPoints, i_IntPoint ;
int iTime, iDof, jDof, iHar, jHar, nVal1, nVal2 = 0, iVal1, iVal2, Type1;
double **H, ***HH, Factor, plus, plus0, weightIntPoint;
int NbrPointsX, OffSet;
struct Expression * Expression_P;
struct Dof * Dofi, *Dofj;
struct Value t_Value;
gMatrix * Jac;
double one=1.0 ;
int iPul, ZeroHarmonic, DcHarmonic;
double E_MH[NBR_MAX_BASISFUNCTIONS][NBR_MAX_BASISFUNCTIONS][NBR_MAX_HARMONIC][NBR_MAX_HARMONIC];
double E_D[NBR_MAX_HARMONIC][NBR_MAX_HARMONIC][MAX_DIM];
void (*xFunctionBFDof[NBR_MAX_BASISFUNCTIONS])
(struct Element * Element, int NumEntity,
double u, double v, double w, double Value[] ) ;
double (*Get_Jacobian)(struct Element*, MATRIX3x3*) ;
void (*Get_IntPoint)(int,int,double*,double*,double*,double*);
double (*Get_Product)(double*,double*,double*) = 0;
/* static double eps; */
GetDP_Begin("Cal_GalerkinTermOfFemEquation_MHJacNL");
FI = EquationTerm_P->Case.LocalTerm.Active ;
QuantityStorage_P = FI->QuantityStorageDof_P ;
/* ---------------------------------------------------------------------- */
/* G e t F u n c t i o n V a l u e f o r t e s t f u n c t i o n s */
/* ---------------------------------------------------------------------- */
if (!(Nbr_Dof = QuantityStorage_P->NbrElementaryBasisFunction)) GetDP_End ;
Get_FunctionValue(Nbr_Dof, (void (**)())xFunctionBFDof,
EquationTerm_P->Case.LocalTerm.Term.TypeOperatorDof,
QuantityStorage_P, &FI->Type_FormDof) ;
for (j = 0 ; j < Nbr_Dof ; j++)
for (k = 0 ; k < Current.NbrHar ; k+=2)
Dof_GetComplexDofValue
(QuantityStorage_P->FunctionSpace->DofData,
QuantityStorage_P->BasisFunction[j].Dof + k/2*gCOMPLEX_INCREMENT,
&Val_Dof[j][k], &Val_Dof[j][k+1]) ;
/* printf("Type1 = %d\n", FI->Type_FormDof); */
nVal1 = NbrValues_Type (Type1 = Get_ValueFromForm(FI->Type_FormDof)) ;
/* ------------------------------------------------------------------- */
/* G e t J a c o b i a n M e t h o d */
/* ------------------------------------------------------------------- */
i = 0 ; while ((i < FI->NbrJacobianCase) &&
((j = (FI->JacobianCase_P0 + i)->RegionIndex) >= 0) &&
!List_Search
(((struct Group *)List_Pointer(Problem_S.Group, j))
->InitialList, &Element->Region, fcmp_int) ) i++ ;
if (i == FI->NbrJacobianCase)
Msg(GERROR, "Undefined Jacobian in Region %d", Element->Region);
Element->JacobianCase = FI->JacobianCase_P0 + i ;
Get_Jacobian = (double (*)(struct Element*, MATRIX3x3*))
Get_JacobianFunction(Element->JacobianCase->TypeJacobian,
Element->Type, &Type_Dimension) ;
if (FI->Flag_ChangeCoord) Get_NodesCoordinatesOfElement(Element) ;
/* integration in space */
IntegrationCase_P =
Get_IntegrationCase(Element, FI->IntegrationCase_L, FI->CriterionIndex);
switch (IntegrationCase_P->Type) {
case ANALYTIC :
Msg(GERROR, "Analytical integration not implemented for MHJacNL");
}
Quadrature_P = (struct Quadrature*)
List_PQuery(IntegrationCase_P->Case, &Element->Type, fcmp_int);
if(!Quadrature_P)
Msg(GERROR, "Unknown type of Element (%s) for Integration method (%s)",
Get_StringForDefine(Element_Type, Element->Type),
((struct IntegrationMethod *)
List_Pointer(Problem_S.IntegrationMethod,
EquationTerm_P->Case.LocalTerm.IntegrationMethodIndex))->Name);
Nbr_IntPoints = Quadrature_P->NumberOfPoints ;
Get_IntPoint = (void(*)(int,int,double*,double*,double*,double*))
Quadrature_P->Function ;
/* integration in fundamental time period */
NbrPointsX = FI->MHJacNL_NbrPointsX;
HH = FI->MHJacNL_HH;
H = FI->MHJacNL_H ;
Expression_P = Problem_Expression0 + FI->MHJacNL_Index;
OffSet = FI->MHJacNL_HarOffSet;
Factor = FI->MHJacNL_Factor;
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
/* C o m p u t a t i o n o f E l e m e n t a r y m a t r i x */
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */
NbrHar = Current.NbrHar ;
/* special treatment of DC-term and associated dummy sinus-term */
DcHarmonic = NbrHar;
ZeroHarmonic = 0;
for (iPul = 0 ; iPul < NbrHar/2 ; iPul++)
if (!Current.DofData->Val_Pulsation[iPul]){
DcHarmonic = 2*iPul ;
ZeroHarmonic = 2*iPul+1 ;
break;
}
/* resetting elementary matrix */
for (iDof = 0 ; iDof < Nbr_Dof ; iDof++)
for (jDof = 0 ; jDof <= iDof ; jDof++)
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = OFFSET ; jHar <= iHar ; jHar++)
E_MH[iDof][jDof][iHar][jHar] = 0. ;
/* volume integration over element */
for (i_IntPoint = 0 ; i_IntPoint < Nbr_IntPoints ; i_IntPoint++) {
Get_IntPoint(Nbr_IntPoints, i_IntPoint,
&Current.u, &Current.v, &Current.w, &weightIntPoint) ;
if (FI->Flag_ChangeCoord) {
Get_BFGeoElement(Element, Current.u, Current.v, Current.w) ;
Element->DetJac = Get_Jacobian(Element, &Element->Jac) ;
weightIntPoint *= fabs(Element->DetJac) ;
if (FI->Flag_InvJac)
Get_InverseMatrix(Type_Dimension, Element->Type, Element->DetJac,
&Element->Jac, &Element->InvJac) ;
}
/* Test and shape Functions (are the same) */
for (i = 0 ; i < Nbr_Dof ; i++) {
xFunctionBFDof[i]
(Element,
QuantityStorage_P->BasisFunction[i].NumEntityInElement+1,
Current.u, Current.v, Current.w, vBFuDof[i]) ;
((void (*)(struct Element*, double*, double*))
FI->xChangeOfCoordinatesEqu) (Element, vBFuDof[i], vBFxDof[i]) ;
}
switch (Type1) {
case SCALAR :
for (k = 0 ; k < NbrHar ; k++){
Val[k] = 0.;
for (j = 0 ; j < Nbr_Dof ; j++)
Val[k] += vBFxDof[j][0] * Val_Dof[j][k] ;
}
break ;
case VECTOR :
for (k = 0 ; k < NbrHar ; k++){
Val[3*k] = Val[3*k+1] = Val[3*k+2] = 0.;
for (j = 0 ; j < Nbr_Dof ; j++){
Val[3*k ] += vBFxDof[j][0] * Val_Dof[j][k] ;
Val[3*k+1] += vBFxDof[j][1] * Val_Dof[j][k] ;
Val[3*k+2] += vBFxDof[j][2] * Val_Dof[j][k] ;
}
}
break ;
}
Current.NbrHar = 1; /* evaluation in time domain */
/* time integration over fundamental period */
for (iTime = 0 ; iTime < NbrPointsX ; iTime++) {
t_Value.Type = Type1 ;
for (iVal1 = 0 ; iVal1 < nVal1 ; iVal1++){
t_Value.Val[iVal1] = 0;
for (iHar = 0 ; iHar < NbrHar ; iHar++)
t_Value.Val[iVal1] += H[iTime][iHar] * Val[iHar*nVal1+iVal1] ;
}
Get_ValueOfExpression(Expression_P, QuantityStorage_P0,
Current.u, Current.v, Current.w, &t_Value);
if (!iTime){
if (!i_IntPoint){
nVal2 = NbrValues_Type (t_Value.Type) ;
Get_Product = (double(*)(double*,double*,double*))
Get_RealProductFunction_Type1xType2xType1 (Type1, t_Value.Type) ;
}
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = OFFSET ; jHar <= iHar ; jHar++)
for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
E_D[iHar][jHar][iVal2] = 0. ;
}
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = OFFSET ; jHar <= iHar ; jHar++){
for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
E_D[iHar][jHar][iVal2] += HH[iTime][iHar][jHar] * t_Value.Val[iVal2] ;
}
} /* for iTime ... */
/*
if (!eps) {
printf("enter value for eps\n");
scanf("%lf",&eps);
printf("eps = %f\n",eps);
}
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = OFFSET ; jHar <= iHar ; jHar++){
for (iVal2 = 0 ; iVal2 < nVal2 ; iVal2++)
if ( E_D[iHar][jHar][iVal2] * E_D[iHar][jHar][iVal2] <
eps * eps * fabs(E_D[iHar][iHar][iVal2] * E_D[jHar][jHar][iVal2]) )
E_D[iHar][jHar][iVal2]=0 ;
}
*/
for (iDof = 0 ; iDof < Nbr_Dof ; iDof++)
for (jDof = 0 ; jDof <= iDof ; jDof++)
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = OFFSET ; jHar <= iHar ; jHar++){
E_MH[iDof][jDof][iHar][jHar] += weightIntPoint *
Get_Product(vBFxDof[iDof], E_D[iHar][jHar], vBFxDof[jDof]) ;
/* printf("%d %d %d %d %e\n", iDof, jDof, iHar, jHar,
E_MH[iDof][jDof][iHar][jHar]) ; */
}
Current.NbrHar = NbrHar ;
} /* for i_IntPoint ... */
/* -------------------------------------------------------------------- */
/* A d d c o n t r i b u t i o n t o J a c o b i a n M a t r i x */
/* -------------------------------------------------------------------- */
Jac = &Current.DofData->Jac;
for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
Dofi = QuantityStorage_P->BasisFunction[iDof].Dof ;
for (jDof = 0 ; jDof <= iDof ; jDof++){
Dofj = QuantityStorage_P->BasisFunction[jDof].Dof ;
for (iHar = 0 ; iHar < NbrHar ; iHar++)
for (jHar = OFFSET ; jHar <= iHar ; jHar++){
plus = plus0 = Factor * E_MH [iDof][jDof] [iHar][jHar] ;
if(jHar==DcHarmonic && iHar!=DcHarmonic) { plus0 *= 1. ; plus *= 2. ;}
Dof_AssembleInMat(Dofi+iHar, Dofj+jHar, 1, &plus, Jac, NULL) ;
if(iHar != jHar)
Dof_AssembleInMat(Dofi+jHar, Dofj+iHar, 1, &plus0, Jac, NULL) ;
if(iDof != jDof){
Dof_AssembleInMat(Dofj+iHar, Dofi+jHar, 1, &plus, Jac, NULL) ;
if(iHar != jHar)
Dof_AssembleInMat(Dofj+jHar, Dofi+iHar, 1, &plus0, Jac, NULL) ;
}
}
}
}
/* dummy 1's on the diagonal for sinus-term of dc-component */
if (ZeroHarmonic) {
for (iDof = 0 ; iDof < Nbr_Dof ; iDof++){
Dofi = QuantityStorage_P->BasisFunction[iDof].Dof + ZeroHarmonic ;
Dof_AssembleInMat(Dofi, Dofi, 1, &one, Jac, NULL) ;
}
}
GetDP_End ;
}
#undef OFFSET
syntax highlighted by Code2HTML, v. 0.9.1