"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "Functions/F_Misc.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.

[ To the main getdp source changes report ]

F_Misc.cpp  (getdp-3.4.0-source.tgz):F_Misc.cpp  (getdp-3.5.0-source.tgz)
// GetDP - Copyright (C) 1997-2021 P. Dular and C. Geuzaine, University of Liege // GetDP - Copyright (C) 1997-2022 P. Dular and C. Geuzaine, University of Liege
// //
// See the LICENSE.txt file for license information. Please report all // See the LICENSE.txt file for license information. Please report all
// issues on https://gitlab.onelab.info/getdp/getdp/issues. // issues on https://gitlab.onelab.info/getdp/getdp/issues.
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <math.h> #include <math.h>
#include "GetDPConfig.h" #include "GetDPConfig.h"
#include "ProData.h" #include "ProData.h"
#include "ProDefine.h" #include "ProDefine.h"
#include "ProParser.h" #include "ProParser.h"
#include "F.h" #include "F.h"
#include "Message.h" #include "Message.h"
#include "OS.h" #include "OS.h"
#include "Cal_Value.h" #include "Cal_Value.h"
#include "Cal_Quantity.h" #include "Cal_Quantity.h"
extern struct CurrentData Current ; extern struct CurrentData Current;
void F_Printf(F_ARG) void F_Printf(F_ARG)
{ {
Print_Value(A) ; Print_Value(A);
printf("\n") ; printf("\n");
} }
void F_Rand(F_ARG) void F_Rand(F_ARG)
{ {
int k; int k;
if(A->Type != SCALAR) if(A->Type != SCALAR)
Message::Error("Non scalar argument for function 'Rand"); Message::Error("Non scalar argument for function 'Rand");
V->Val[0] = A->Val[0] * (double)rand() / (double)RAND_MAX; V->Val[0] = A->Val[0] * (double)rand() / (double)RAND_MAX;
if (Current.NbrHar != 1){ if(Current.NbrHar != 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR ; V->Type = SCALAR;
} }
void F_CompElementNum (F_ARG) void F_CompElementNum(F_ARG)
{ {
if(!Current.Element || !Current.ElementSource) if(!Current.Element || !Current.ElementSource)
Message::Error("Uninitialized Element in 'F_CompElementNum'"); Message::Error("Uninitialized Element in 'F_CompElementNum'");
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = (Current.Element->Num == Current.ElementSource->Num) ; V->Val[0] = (Current.Element->Num == Current.ElementSource->Num);
} }
void F_ElementNum (F_ARG) void F_ElementNum(F_ARG)
{ {
if(!Current.Element) if(!Current.Element)
Message::Error("Uninitialized Element in 'F_ElementNum'"); Message::Error("Uninitialized Element in 'F_ElementNum'");
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = Current.Element->Num ; V->Val[0] = Current.Element->Num;
} }
void F_QuadraturePointIndex (F_ARG) void F_QuadraturePointIndex(F_ARG)
{ {
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = Current.QuadraturePointIndex ; V->Val[0] = Current.QuadraturePointIndex;
} }
void F_GetCpuTime (F_ARG) void F_GetCpuTime(F_ARG)
{ {
double s = 0.; double s = 0.;
long mem = 0; long mem = 0;
GetResources(&s, &mem); GetResources(&s, &mem);
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = s ; V->Val[0] = s;
} }
void F_GetWallClockTime (F_ARG) void F_GetWallClockTime(F_ARG)
{ {
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = Message::GetWallClockTime() ; V->Val[0] = Message::GetWallClockTime();
} }
void F_GetMemory (F_ARG) void F_GetMemory(F_ARG)
{ {
double s = 0.; double s = 0.;
long mem = 0; long mem = 0;
GetResources(&s, &mem); GetResources(&s, &mem);
double val = mem / 1024. / 1024.; double val = mem / 1024. / 1024.;
V->Type = SCALAR ; V->Type = SCALAR;
V->Val[0] = val ; V->Val[0] = val;
} }
void F_SetNumberRunTime (F_ARG) void F_SetNumberRunTime(F_ARG)
{ {
double val = A->Val[0]; double val = A->Val[0];
int type = A->Type; int type = A->Type;
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++) V->Val[MAX_DIM * k] = 0.;
V->Val[MAX_DIM * k] = 0. ;
V->Type = SCALAR; V->Type = SCALAR;
if(type != SCALAR){ if(type != SCALAR) {
Message::Error("Non scalar argument for function 'SetNumberRunTime'"); Message::Error("Non scalar argument for function 'SetNumberRunTime'");
return; return;
} }
if(!Fct->String){ if(!Fct->String) {
Message::Error("Missing ONELAB variable name: use SetNumberRunTime[value]{\" Message::Error(
name\"}"); "Missing ONELAB variable name: use SetNumberRunTime[value]{\"name\"}");
return; return;
} }
Message::SetOnelabNumber(Fct->String, val); Message::SetOnelabNumber(Fct->String, val);
V->Val[0] = val ; V->Val[0] = val;
} }
void F_SetNumberRunTimeWithChoices (F_ARG) void F_SetNumberRunTimeWithChoices(F_ARG)
{ {
double val = A->Val[0]; double val = A->Val[0];
int type = A->Type; int type = A->Type;
for (int k = 0; k < Current.NbrHar; k++) for(int k = 0; k < Current.NbrHar; k++) V->Val[MAX_DIM * k] = 0.;
V->Val[MAX_DIM * k] = 0. ;
V->Type = SCALAR; V->Type = SCALAR;
if(type != SCALAR){ if(type != SCALAR) {
Message::Error("Non scalar argument for function 'SetNumberRunTime'"); Message::Error("Non scalar argument for function 'SetNumberRunTime'");
return; return;
} }
if(!Fct->String){ if(!Fct->String) {
Message::Error("Missing ONELAB variable name: use SetNumberRunTime[value]{\" Message::Error(
name\"}"); "Missing ONELAB variable name: use SetNumberRunTime[value]{\"name\"}");
return; return;
} }
std::vector<double> v(1, val); std::vector<double> v(1, val);
Message::AddOnelabNumberChoice(Fct->String, v); Message::AddOnelabNumberChoice(Fct->String, v);
V->Val[0] = val ; V->Val[0] = val;
} }
void F_GetNumberRunTime (F_ARG) void F_GetNumberRunTime(F_ARG)
{ {
if(Fct->NbrArguments){ if(Fct->NbrArguments) { Cal_CopyValue(A, V); }
Cal_CopyValue(A, V); else {
} for(int k = 0; k < Current.NbrHar; k++) V->Val[MAX_DIM * k] = 0.;
else{
for (int k = 0; k < Current.NbrHar; k++)
V->Val[MAX_DIM * k] = 0. ;
V->Type = SCALAR; V->Type = SCALAR;
} }
if(!Fct->String){ if(!Fct->String) {
Message::Error("Missing ONELAB variable name: use GetNumberRunTime[]{\"name\ Message::Error(
"}"); "Missing ONELAB variable name: use GetNumberRunTime[]{\"name\"}");
return; return;
} }
if(Message::UseOnelab()) if(Message::UseOnelab()) V->Val[0] = Message::GetOnelabNumber(Fct->String);
V->Val[0] = Message::GetOnelabNumber(Fct->String);
} }
void F_SetVariable (F_ARG) void F_SetVariable(F_ARG)
{ {
if(!Fct->String){ if(!Fct->String) {
Message::Error("Missing runtime variable name: use SetVariable[...]{$name}") Message::Error(
; "Missing runtime variable name: use SetVariable[...]{$name}");
return; return;
} }
if(Fct->NbrArguments < 1){ if(Fct->NbrArguments < 1) {
Message::Error("Missing argument in SetVariable[...]{$name}"); Message::Error("Missing argument in SetVariable[...]{$name}");
return; return;
} }
Cal_CopyValue(A, V); Cal_CopyValue(A, V);
char tmp[256]; char tmp[256];
strcpy(tmp, Fct->String); strcpy(tmp, Fct->String);
for(int i = 1; i < Fct->NbrArguments; i++){ for(int i = 1; i < Fct->NbrArguments; i++) {
if((A+i)->Type != SCALAR){ if((A + i)->Type != SCALAR) {
Message::Error("Non-scalar argument in SetVariable"); Message::Error("Non-scalar argument in SetVariable");
return; return;
} }
char tmp2[256]; char tmp2[256];
sprintf(tmp2, "_%g", (A+i)->Val[0]); sprintf(tmp2, "_%g", (A + i)->Val[0]);
strcat(tmp, tmp2); strcat(tmp, tmp2);
} }
Cal_StoreInVariable(V, tmp); Cal_StoreInVariable(V, tmp);
} }
void F_SetCumulativeVariable (F_ARG) void F_SetCumulativeVariable(F_ARG)
{ {
if(Fct->NbrArguments < 1){ if(Fct->NbrArguments < 1) {
Message::Error("Missing argument in SetCumulativeVariable[...]{$name}"); Message::Error("Missing argument in SetCumulativeVariable[...]{$name}");
return; return;
} }
Cal_CopyValue(A, V); Cal_CopyValue(A, V);
char tmp[256]; char tmp[256];
strcpy(tmp, Fct->String); strcpy(tmp, Fct->String);
for(int i = 1; i < Fct->NbrArguments; i++){ for(int i = 1; i < Fct->NbrArguments; i++) {
if((A+i)->Type != SCALAR){ if((A + i)->Type != SCALAR) {
Message::Error("Non-scalar argument in SetCumulativeVariable"); Message::Error("Non-scalar argument in SetCumulativeVariable");
return; return;
} }
char tmp2[256]; char tmp2[256];
sprintf(tmp2, "_%g", (A+i)->Val[0]); sprintf(tmp2, "_%g", (A + i)->Val[0]);
strcat(tmp, tmp2); strcat(tmp, tmp2);
} }
struct Value V_saved ; struct Value V_saved;
Cal_GetValueSaved(&V_saved, tmp); Cal_GetValueSaved(&V_saved, tmp);
Cal_AddValue(V, &V_saved, V); Cal_AddValue(V, &V_saved, V);
Cal_StoreInVariable(V, tmp); Cal_StoreInVariable(V, tmp);
} }
void F_GetVariable (F_ARG) void F_GetVariable(F_ARG)
{ {
if(!Fct->String){ if(!Fct->String) {
Message::Error("Missing runtime variable name: use GetVariable[...]{$name}") Message::Error(
; "Missing runtime variable name: use GetVariable[...]{$name}");
return; return;
} }
char tmp[256]; char tmp[256];
strcpy(tmp, Fct->String); strcpy(tmp, Fct->String);
for(int i = 0; i < Fct->NbrArguments; i++){ for(int i = 0; i < Fct->NbrArguments; i++) {
if((A+i)->Type != SCALAR){ if((A + i)->Type != SCALAR) {
Message::Error("Non-scalar argument in GetVariable"); Message::Error("Non-scalar argument in GetVariable");
return; return;
} }
char tmp2[256]; char tmp2[256];
sprintf(tmp2, "_%g", (A+i)->Val[0]); sprintf(tmp2, "_%g", (A + i)->Val[0]);
strcat(tmp, tmp2); strcat(tmp, tmp2);
} }
Cal_GetValueSaved(V, tmp); Cal_GetValueSaved(V, tmp);
} }
void F_ValueFromTable (F_ARG) void F_ValueFromTable(F_ARG)
{ {
if(!Fct->String){ if(!Fct->String) {
Message::Error("Missing ElementTable or NodeTable name: use ValueFromTable[. Message::Error(
..]{name}"); "Missing ElementTable or NodeTable name: use ValueFromTable[...]{name}");
return; return;
} }
std::map<int, std::vector<double> > &table(GetDPNumbersMap[Fct->String]); std::map<int, std::vector<double> > &table(GetDPNumbersMap[Fct->String]);
std::vector<double> &val(table[Current.NumEntity]); std::vector<double> &val(table[Current.NumEntity]);
if(val.size() == 1){ if(val.size() == 1) {
V->Val[0] = val[0]; V->Val[0] = val[0];
V->Type = SCALAR ; V->Type = SCALAR;
return; return;
} }
//if(GetDPNumbersMap.size()){ // if(GetDPNumbersMap.size()){
// Message::Warning("No element or node table found with name %s, or " // Message::Warning("No element or node table found with name %s, or "
// "no entity index %d in the table", Fct->String, // "no entity index %d in the table", Fct->String,
// Current.NumEntity); // Current.NumEntity);
//} // }
for(int i = 0; i < Fct->NbrArguments; i++){ for(int i = 0; i < Fct->NbrArguments; i++) {
if(i != 0){ if(i != 0) {
Message::Warning("Ignoring %d-th argument in ValueFromTable"); Message::Warning("Ignoring %d-th argument in ValueFromTable");
continue; continue;
} }
if((A+i)->Type != SCALAR){ if((A + i)->Type != SCALAR) {
Message::Error("Non-scalar default argument in ValueFromTable"); Message::Error("Non-scalar default argument in ValueFromTable");
return; return;
} }
Cal_CopyValue(A + i, V); Cal_CopyValue(A + i, V);
return; return;
} }
Message::Warning("Missing table data or default value in ValueFromTable"); Message::Warning("Missing table data or default value in ValueFromTable");
V->Val[0] = 0. ; V->Val[0] = 0.;
V->Type = SCALAR ; V->Type = SCALAR;
} }
void F_VirtualWork (F_ARG) void F_VirtualWork(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element)
Message::Error("Uninitialized Element in 'F_VirtualWork'"); Message::Error("Uninitialized Element in 'F_VirtualWork'");
Current.flagAssDiag = 1; /*+++prov*/ Current.flagAssDiag = 1; /*+++prov*/
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) { if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
s[0] = ( DetJac_dx[0] * ( - squF[0] + squF[1] + squF[2] ) s[0] = (DetJac_dx[0] * (-squF[0] + squF[1] + squF[2]) -
- 2 * DetJac_dx[1] * squF[3] 2 * DetJac_dx[1] * squF[3] - 2 * DetJac_dx[2] * squF[5]) /
- 2 * DetJac_dx[2] * squF[5])/DetJac ; DetJac;
s[1] = ( DetJac_dx[1] * ( squF[0] - squF[1] + squF[2] ) s[1] = (DetJac_dx[1] * (squF[0] - squF[1] + squF[2]) -
- 2 * DetJac_dx[0] * squF[3] 2 * DetJac_dx[0] * squF[3] - 2 * DetJac_dx[2] * squF[4]) /
- 2 * DetJac_dx[2] * squF[4])/DetJac ; DetJac;
s[2] = ( DetJac_dx[2] * ( squF[0] + squF[1] - squF[2] ) s[2] = (DetJac_dx[2] * (squF[0] + squF[1] - squF[2]) -
- 2 * DetJac_dx[0] * squF[5] 2 * DetJac_dx[0] * squF[5] - 2 * DetJac_dx[1] * squF[4]) /
- 2 * DetJac_dx[1] * squF[4])/DetJac ; DetJac;
} }
else { else {
Message::Warning("Zero determinant in 'F_VirtualWork'") ; Message::Warning("Zero determinant in 'F_VirtualWork'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_NodeForceDensity (F_ARG) void F_NodeForceDensity(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double Grad_n[3] ; double Grad_n[3];
double s11 = 0., s12 = 0., s13 = 0. ; double s11 = 0., s12 = 0., s13 = 0.;
double s21 = 0., s22 = 0., s23 = 0. ; double s21 = 0., s22 = 0., s23 = 0.;
double s31 = 0., s32 = 0., s33 = 0. ; double s31 = 0., s32 = 0., s33 = 0.;
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element)
Message::Error("Uninitialized Element in 'F_NodeForceDensity'"); Message::Error("Uninitialized Element in 'F_NodeForceDensity'");
Current.flagAssDiag = 1; /*+++prov*/ Current.flagAssDiag = 1; /*+++prov*/
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(A->Type == TENSOR_SYM){ if(i < Current.Element->GeoElement->NbrNodes) {
s11 = A->Val[0] ; if(A->Type == TENSOR_SYM) {
s12 = A->Val[1] ; s11 = A->Val[0];
s13 = A->Val[2] ; s12 = A->Val[1];
s13 = A->Val[2];
s21 = s12; s21 = s12;
s22 = A->Val[3] ; s22 = A->Val[3];
s23 = A->Val[4] ; s23 = A->Val[4];
s31 = s13; s31 = s13;
s32 = s23; s32 = s23;
s33 = A->Val[5] ; s33 = A->Val[5];
} }
else if(A->Type == TENSOR){ else if(A->Type == TENSOR) {
s11 = A->Val[0] ; s11 = A->Val[0];
s12 = A->Val[1] ; s12 = A->Val[1];
s13 = A->Val[2] ; s13 = A->Val[2];
s21 = A->Val[3] ; s21 = A->Val[3];
s22 = A->Val[4] ; s22 = A->Val[4];
s23 = A->Val[5] ; s23 = A->Val[5];
s31 = A->Val[6] ; s31 = A->Val[6];
s32 = A->Val[7] ; s32 = A->Val[7];
s33 = A->Val[8] ; s33 = A->Val[8];
} }
else{ else {
Message::Error("Unknown tensor type in 'F_NodeForceDensity'") ; Message::Error("Unknown tensor type in 'F_NodeForceDensity'");
} }
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
Grad_n[0] = Grad_n[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
Grad_n[1] = Grad_n[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
Grad_n[2] = Grad_n[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
s[0] = ( Grad_n[0] * s11 + Grad_n[1] * s12 + Grad_n[2] * s13 ) / DetJac ; s[0] = (Grad_n[0] * s11 + Grad_n[1] * s12 + Grad_n[2] * s13) / DetJac;
s[1] = ( Grad_n[0] * s21 + Grad_n[1] * s22 + Grad_n[2] * s23 ) / DetJac ; s[1] = (Grad_n[0] * s21 + Grad_n[1] * s22 + Grad_n[2] * s23) / DetJac;
s[2] = ( Grad_n[0] * s31 + Grad_n[1] * s32 + Grad_n[2] * s33 ) / DetJac ; s[2] = (Grad_n[0] * s31 + Grad_n[1] * s32 + Grad_n[2] * s33) / DetJac;
} }
else { else {
Message::Warning("Zero determinant in 'F_NodeForceDensity'") ; Message::Warning("Zero determinant in 'F_NodeForceDensity'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
// Blex added 25/04/14 update 06/06/14 // Blex added 25/04/14 update 06/06/14
void F_Felec (F_ARG) void F_Felec(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_Felec'");
Message::Error("Uninitialized Element in 'F_Felec'");
Current.flagAssDiag = 1; /*+++prov*/ Current.flagAssDiag = 1; /*+++prov*/
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - 0.5 * ( dsdx[0] * ( squF[0] - squF[1] - squF[2] ) s[0] = -0.5 * (dsdx[0] * (squF[0] - squF[1] - squF[2]) +
+ 2 * dsdx[1] * squF[3] 2 * dsdx[1] * squF[3] + 2 * dsdx[2] * squF[5]);
+ 2 * dsdx[2] * squF[5]) ;
s[1] = -0.5 * (dsdx[1] * (-squF[0] + squF[1] - squF[2]) +
s[1] = - 0.5 * ( dsdx[1] * ( - squF[0] + squF[1] - squF[2] ) 2 * dsdx[0] * squF[3] + 2 * dsdx[2] * squF[4]);
+ 2 * dsdx[0] * squF[3]
+ 2 * dsdx[2] * squF[4]) ; s[2] = -0.5 * (dsdx[2] * (-squF[0] - squF[1] + squF[2]) +
2 * dsdx[0] * squF[5] + 2 * dsdx[1] * squF[4]);
s[2] = - 0.5 * ( dsdx[2] * ( - squF[0] - squF[1] + squF[2] )
+ 2 * dsdx[0] * squF[5]
+ 2 * dsdx[1] * squF[4]) ;
} }
else { else {
Message::Warning("Zero determinant in 'F_Felec'") ; Message::Warning("Zero determinant in 'F_Felec'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFxdux (F_ARG) void F_dFxdux(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFxdux'");
Message::Error("Uninitialized Element in 'F_dFxdux'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - squF[0] * dsdx[0] ; s[0] = -squF[0] * dsdx[0];
s[1] = - squF[0] * dsdx[1] ; s[1] = -squF[0] * dsdx[1];
s[2] = - squF[0] * dsdx[2] ; s[2] = -squF[0] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dFxdux'") ; Message::Warning("Zero determinant in 'F_dFxdux'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFydux (F_ARG) void F_dFydux(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFydux'");
Message::Error("Uninitialized Element in 'F_dFydux'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = -0.5 *
(2 * squF[3] * dsdx[0] + (-squF[0] - squF[1] + squF[2]) * dsdx[1] -
2 * squF[4] * dsdx[2]);
s[0] = - 0.5 * (2 * squF[3] * dsdx[0] + (- squF[0] - squF[1] + squF[2]) * s[1] = -0.5 * ((squF[0] + squF[1] - squF[2]) * dsdx[0] +
dsdx[1] - 2 * squF[4] * dsdx[2]) ; 2 * squF[3] * dsdx[1] + 2 * squF[5] * dsdx[2]);
s[1] = - 0.5 * ((squF[0] + squF[1] - squF[2]) * dsdx[0] + 2 * squF[3] * ds s[2] = -0.5 * (2 * squF[4] * dsdx[0] - 2 * squF[5] * dsdx[1] +
dx[1] + 2 * squF[5] * dsdx[2]) ; 2 * squF[3] * dsdx[2]);
s[2] = - 0.5 * (2 * squF[4] * dsdx[0] - 2 * squF[5] * dsdx[1] + 2 * squF[3
] * dsdx[2]) ;
} }
else { else {
Message::Warning("Zero determinant in 'F_dFydux'") ; Message::Warning("Zero determinant in 'F_dFydux'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFzdux (F_ARG) void F_dFzdux(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFzdux'");
Message::Error("Uninitialized Element in 'F_dFzdux'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0) {
dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2] / DetJac;
if(DetJac != 0){ s[0] = -0.5 * (2 * squF[5] * dsdx[0] - 2 * squF[4] * dsdx[1] +
dsdx[0] = DetJac_dx[0]/DetJac ; (-squF[0] + squF[1] - squF[2]) * dsdx[2]);
dsdx[1] = DetJac_dx[1]/DetJac ;
dsdx[2] = DetJac_dx[2]/DetJac ;
s[0] = - 0.5 * (2 * squF[5] * dsdx[0] - 2 * squF[4] * dsdx[1] + (-squF[0] s[1] = -0.5 * (2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] -
+ squF[1] - squF[2]) * dsdx[2]) ; 2 * squF[3] * dsdx[2]);
s[1] = - 0.5 * (2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] - 2 * squF[3 s[2] = -0.5 * ((squF[0] - squF[1] + squF[2]) * dsdx[0] +
] * dsdx[2]) ; 2 * squF[3] * dsdx[1] + 2 * squF[5] * dsdx[2]);
s[2] = - 0.5 * ((squF[0] - squF[1] + squF[2]) * dsdx[0] + 2 * squF[3] * ds
dx[1] + 2 * squF[5] * dsdx[2]) ;
} }
else { else {
Message::Warning("Zero determinant in 'F_dFzdux'") ; Message::Warning("Zero determinant in 'F_dFzdux'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFxduy (F_ARG) void F_dFxduy(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFxduy'");
Message::Error("Uninitialized Element in 'F_dFxduy'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0) {
dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = -0.5 *
(2 * squF[3] * dsdx[0] + (squF[0] + squF[1] - squF[2]) * dsdx[1] +
2 * squF[4] * dsdx[2]);
if(DetJac != 0){ s[1] = -0.5 * ((-squF[0] - squF[1] + squF[2]) * dsdx[0] +
dsdx[0] = DetJac_dx[0]/DetJac ; 2 * squF[3] * dsdx[1] - 2 * squF[5] * dsdx[2]);
dsdx[1] = DetJac_dx[1]/DetJac ;
dsdx[2] = DetJac_dx[2]/DetJac ;
s[0] = - 0.5 * (2 * squF[3] * dsdx[0] + (squF[0] + squF[1] - squF[2]) * ds s[2] = -0.5 * (-2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] +
dx[1] + 2 * squF[4] * dsdx[2]) ; 2 * squF[3] * dsdx[2]);
s[1] = - 0.5 * ((-squF[0] - squF[1] + squF[2]) * dsdx[0] + 2 * squF[3] * d
sdx[1] - 2 * squF[5] * dsdx[2]) ;
s[2] = - 0.5 * (- 2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] + 2 * squF
[3] * dsdx[2]) ;
} }
else { else {
Message::Warning("Zero determinant in 'F_dFxduy'") ; Message::Warning("Zero determinant in 'F_dFxduy'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFyduy (F_ARG) void F_dFyduy(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFyduy'");
Message::Error("Uninitialized Element in 'F_dFyduy'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - squF[1] * dsdx[0] ; s[0] = -squF[1] * dsdx[0];
s[1] = - squF[1] * dsdx[1] ; s[1] = -squF[1] * dsdx[1];
s[2] = - squF[1] * dsdx[2] ; s[2] = -squF[1] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dFyduy'") ; Message::Warning("Zero determinant in 'F_dFyduy'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFzduy (F_ARG) void F_dFzduy(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFzduy'");
Message::Error("Uninitialized Element in 'F_dFzduy'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - 0.5 * (2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] - 2 * squF[3 s[0] = -0.5 * (2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] -
] * dsdx[2]) ; 2 * squF[3] * dsdx[2]);
s[1] = - 0.5 * (-2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] + (squF[0] s[1] = -0.5 * (-2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] +
- squF[1] - squF[2]) * dsdx[2]) ; (squF[0] - squF[1] - squF[2]) * dsdx[2]);
s[2] = - 0.5 * (2 * squF[3] * dsdx[0] + (-squF[0] + squF[1] + squF[2]) * d s[2] = -0.5 *
sdx[1] + 2 * squF[4] * dsdx[2]) ; (2 * squF[3] * dsdx[0] + (-squF[0] + squF[1] + squF[2]) * dsdx[1] +
2 * squF[4] * dsdx[2]);
} }
else { else {
Message::Warning("Zero determinant in 'F_dFzduy'") ; Message::Warning("Zero determinant in 'F_dFzduy'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFxduz (F_ARG) void F_dFxduz(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFxduz'");
Message::Error("Uninitialized Element in 'F_dFxduz'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - 0.5 * (2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] + (squF[0] - s[0] = -0.5 * (2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] +
squF[1] + squF[2]) * dsdx[2]) ; (squF[0] - squF[1] + squF[2]) * dsdx[2]);
s[1] = - 0.5 * (-2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] + 2 * squF[ s[1] = -0.5 * (-2 * squF[4] * dsdx[0] + 2 * squF[5] * dsdx[1] +
3] * dsdx[2]) ; 2 * squF[3] * dsdx[2]);
s[2] = - 0.5 * ((-squF[0] + squF[1] - squF[2]) * dsdx[0] - 2 * squF[3] * d s[2] = -0.5 * ((-squF[0] + squF[1] - squF[2]) * dsdx[0] -
sdx[1] + 2 * squF[5] * dsdx[2]) ; 2 * squF[3] * dsdx[1] + 2 * squF[5] * dsdx[2]);
} }
else { else {
Message::Warning("Zero determinant in 'F_dFxduz'") ; Message::Warning("Zero determinant in 'F_dFxduz'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFyduz (F_ARG) void F_dFyduz(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFyduz'");
Message::Error("Uninitialized Element in 'F_dFyduz'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - 0.5 * (2 * squF[4] * dsdx[0] - 2 * squF[5] * dsdx[1] + 2 * squF[3 s[0] = -0.5 * (2 * squF[4] * dsdx[0] - 2 * squF[5] * dsdx[1] +
] * dsdx[2]) ; 2 * squF[3] * dsdx[2]);
s[1] = - 0.5 * (2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] + (-squF[0] s[1] = -0.5 * (2 * squF[5] * dsdx[0] + 2 * squF[4] * dsdx[1] +
+ squF[1] + squF[2]) * dsdx[2]) ; (-squF[0] + squF[1] + squF[2]) * dsdx[2]);
s[2] = - 0.5 * (-2 * squF[3] * dsdx[0] + (squF[0] - squF[1] - squF[2]) * d s[2] = -0.5 *
sdx[1] + 2 * squF[4] * dsdx[2]) ; (-2 * squF[3] * dsdx[0] + (squF[0] - squF[1] - squF[2]) * dsdx[1] +
2 * squF[4] * dsdx[2]);
} }
else { else {
Message::Warning("Zero determinant in 'F_dFyduz'") ; Message::Warning("Zero determinant in 'F_dFyduz'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFzduz (F_ARG) void F_dFzduz(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[6] ; double DetJac_dx[3], squF[6];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFzduz'");
Message::Error("Uninitialized Element in 'F_dFzduz'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) {
squF[j] = A->Val[j] * A->Val[j] ; squF[j] = A->Val[j] * A->Val[j];
squF[j+3] = A->Val[j] * A->Val[(j<2)?j+1:0] ; squF[j + 3] = A->Val[j] * A->Val[(j < 2) ? j + 1 : 0];
} }
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac ; DetJac = Current.Element->DetJac;
Jac = Current.Element->Jac ; Jac = Current.Element->Jac;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - squF[2] * dsdx[0] ; s[0] = -squF[2] * dsdx[0];
s[1] = - squF[2] * dsdx[1] ; s[1] = -squF[2] * dsdx[1];
s[2] = - squF[2] * dsdx[2] ; s[2] = -squF[2] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dFzduz'") ; Message::Warning("Zero determinant in 'F_dFzduz'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFxdv (F_ARG) void F_dFxdv(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[3] ; double DetJac_dx[3], squF[3];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFxdv'");
Message::Error("Uninitialized Element in 'F_dFxdv'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) { squF[j] = A->Val[j]; }
squF[j] = A->Val[j] ;
} // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac;
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); Jac = Current.Element->Jac;
DetJac = Current.Element->DetJac ;
Jac = Current.Element->Jac ;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2] ; s[0] = -squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2];
s[1] = squF[1] * dsdx[0] - squF[0] * dsdx[1] ; s[1] = squF[1] * dsdx[0] - squF[0] * dsdx[1];
s[2] = squF[2] * dsdx[0] - squF[0] * dsdx[2] ; s[2] = squF[2] * dsdx[0] - squF[0] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dFxdv'") ; Message::Warning("Zero determinant in 'F_dFxdv'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFydv (F_ARG) void F_dFydv(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[3] ; double DetJac_dx[3], squF[3];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFydv'");
Message::Error("Uninitialized Element in 'F_dFydv'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) { squF[j] = A->Val[j]; }
squF[j] = A->Val[j] ;
} // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac;
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); Jac = Current.Element->Jac;
DetJac = Current.Element->DetJac ;
Jac = Current.Element->Jac ;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - squF[1] * dsdx[0] + squF[0] * dsdx[1] ; s[0] = -squF[1] * dsdx[0] + squF[0] * dsdx[1];
s[1] = - squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2] ; s[1] = -squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2];
s[2] = squF[2] * dsdx[1] - squF[1] * dsdx[2] ; s[2] = squF[2] * dsdx[1] - squF[1] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dFydv'") ; Message::Warning("Zero determinant in 'F_dFydv'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dFzdv (F_ARG) void F_dFzdv(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[3] ; double DetJac_dx[3], squF[3];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dFzdv'");
Message::Error("Uninitialized Element in 'F_dFzdv'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) { squF[j] = A->Val[j]; }
squF[j] = A->Val[j] ;
} // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac;
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); Jac = Current.Element->Jac;
DetJac = Current.Element->DetJac ;
Jac = Current.Element->Jac ;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = - squF[2] * dsdx[0] + squF[0] * dsdx[2] ; s[0] = -squF[2] * dsdx[0] + squF[0] * dsdx[2];
s[1] = - squF[2] * dsdx[1] + squF[1] * dsdx[2] ; s[1] = -squF[2] * dsdx[1] + squF[1] * dsdx[2];
s[2] = - squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2] ; s[2] = -squF[0] * dsdx[0] - squF[1] * dsdx[1] - squF[2] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dFzdv'") ; Message::Warning("Zero determinant in 'F_dFzdv'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dWedxdv (F_ARG) void F_dWedxdv(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[3] ; double DetJac_dx[3], squF[3];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dWedxdv'");
Message::Error("Uninitialized Element in 'F_dWedxdv'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) { squF[j] = A->Val[j]; }
squF[j] = A->Val[j] ;
} // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac;
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); Jac = Current.Element->Jac;
DetJac = Current.Element->DetJac ;
Jac = Current.Element->Jac ;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = -squF[0] * dsdx[0] + squF[1] * dsdx[1] + squF[2] * dsdx[2] ; s[0] = -squF[0] * dsdx[0] + squF[1] * dsdx[1] + squF[2] * dsdx[2];
s[1] = -squF[1] * dsdx[0] - squF[0] * dsdx[1] ; s[1] = -squF[1] * dsdx[0] - squF[0] * dsdx[1];
s[2] = -squF[2] * dsdx[0] - squF[0] * dsdx[2] ; s[2] = -squF[2] * dsdx[0] - squF[0] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dWedxdv'") ; Message::Warning("Zero determinant in 'F_dWedxdv'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dWedydv (F_ARG) void F_dWedydv(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[3] ; double DetJac_dx[3], squF[3];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dWedydv'");
Message::Error("Uninitialized Element in 'F_dWedydv'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) { squF[j] = A->Val[j]; }
squF[j] = A->Val[j] ;
} // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac;
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); Jac = Current.Element->Jac;
DetJac = Current.Element->DetJac ;
Jac = Current.Element->Jac ;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = -squF[1] * dsdx[0] - squF[0] * dsdx[1] ; s[0] = -squF[1] * dsdx[0] - squF[0] * dsdx[1];
s[1] = squF[0] * dsdx[0] - squF[1] * dsdx[1] + squF[2] * dsdx[2] ; s[1] = squF[0] * dsdx[0] - squF[1] * dsdx[1] + squF[2] * dsdx[2];
s[2] = -squF[2] * dsdx[1] - squF[1] * dsdx[2] ; s[2] = -squF[2] * dsdx[1] - squF[1] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dWedydv'") ; Message::Warning("Zero determinant in 'F_dWedydv'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
void F_dWedzdv (F_ARG) void F_dWedzdv(F_ARG)
{ {
MATRIX3x3 Jac ; MATRIX3x3 Jac;
double DetJac ; double DetJac;
double DetJac_dx[3], squF[3] ; double DetJac_dx[3], squF[3];
double dsdx[3] ; //Derivative of the base functions with respect to x, y and z double dsdx[3]; // Derivative of the base functions with respect to x, y and z
double s[3] = {0.,0.,0.}; double s[3] = {0., 0., 0.};
if(!Current.Element) if(!Current.Element) Message::Error("Uninitialized Element in 'F_dWedzdv'");
Message::Error("Uninitialized Element in 'F_dWedzdv'");
int numNode = Current.NumEntity; int numNode = Current.NumEntity;
int i = 0 ; int i = 0;
while (i < Current.Element->GeoElement->NbrNodes && while(i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i] != numNode) i++; Current.Element->GeoElement->NumNodes[i] != numNode)
i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
if(i < Current.Element->GeoElement->NbrNodes) {
for(int j = 0; j < 3; j++) { for(int j = 0; j < 3; j++) { squF[j] = A->Val[j]; }
squF[j] = A->Val[j] ;
} // Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w);
DetJac = Current.Element->DetJac;
//Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w); Jac = Current.Element->Jac;
DetJac = Current.Element->DetJac ;
Jac = Current.Element->Jac ;
DetJac_dx[0] = DetJac_dx[0] =
Current.Element->dndu[i][0] * ( Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32 ) Current.Element->dndu[i][0] * (Jac.c22 * Jac.c33 - Jac.c23 * Jac.c32) -
- Current.Element->dndu[i][1] * ( Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32 ) Current.Element->dndu[i][1] * (Jac.c12 * Jac.c33 - Jac.c13 * Jac.c32) +
+ Current.Element->dndu[i][2] * ( Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13 ); Current.Element->dndu[i][2] * (Jac.c12 * Jac.c23 - Jac.c22 * Jac.c13);
DetJac_dx[1] = DetJac_dx[1] =
- Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31 ) -Current.Element->dndu[i][0] * (Jac.c21 * Jac.c33 - Jac.c23 * Jac.c31) +
+ Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c33 - Jac.c13 * Jac.c31) -
- Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c23 - Jac.c13 * Jac.c21);
DetJac_dx[2] = DetJac_dx[2] =
Current.Element->dndu[i][0] * ( Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31 ) Current.Element->dndu[i][0] * (Jac.c21 * Jac.c32 - Jac.c22 * Jac.c31) -
- Current.Element->dndu[i][1] * ( Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31 ) Current.Element->dndu[i][1] * (Jac.c11 * Jac.c32 - Jac.c12 * Jac.c31) +
+ Current.Element->dndu[i][2] * ( Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21 ); Current.Element->dndu[i][2] * (Jac.c11 * Jac.c22 - Jac.c12 * Jac.c21);
if(DetJac != 0){ if(DetJac != 0) {
dsdx[0] = DetJac_dx[0]/DetJac ; dsdx[0] = DetJac_dx[0] / DetJac;
dsdx[1] = DetJac_dx[1]/DetJac ; dsdx[1] = DetJac_dx[1] / DetJac;
dsdx[2] = DetJac_dx[2]/DetJac ; dsdx[2] = DetJac_dx[2] / DetJac;
s[0] = -squF[2] * dsdx[0] - squF[0] * dsdx[2] ; s[0] = -squF[2] * dsdx[0] - squF[0] * dsdx[2];
s[1] = -squF[2] * dsdx[1] - squF[1] * dsdx[2] ; s[1] = -squF[2] * dsdx[1] - squF[1] * dsdx[2];
s[2] = squF[0] * dsdx[0] + squF[1] * dsdx[1] - squF[2] * dsdx[2] ; s[2] = squF[0] * dsdx[0] + squF[1] * dsdx[1] - squF[2] * dsdx[2];
} }
else { else {
Message::Warning("Zero determinant in 'F_dWedzdv'") ; Message::Warning("Zero determinant in 'F_dWedzdv'");
} }
} }
V->Type = VECTOR ; V->Type = VECTOR;
V->Val[0] = s[0] ; V->Val[0] = s[0];
V->Val[1] = s[1] ; V->Val[1] = s[1];
V->Val[2] = s[2] ; V->Val[2] = s[2];
} }
// End Blex added // End Blex added
void F_AssDiag(F_ARG) void F_AssDiag(F_ARG)
{ {
int k ; int k;
if (Fct->NbrParameters == 1) if(Fct->NbrParameters == 1)
Current.flagAssDiag = Fct->Para[0]; Current.flagAssDiag = Fct->Para[0];
else else
Current.flagAssDiag = 2; /*+++prov*/ Current.flagAssDiag = 2; /*+++prov*/
V->Val[0] = 1.; V->Val[0] = 1.;
if (Current.NbrHar != 1){ if(Current.NbrHar != 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
V->Type = SCALAR ; V->Type = SCALAR;
} }
void F_AtIndex(F_ARG) void F_AtIndex(F_ARG)
{ {
int index = 0; int index = 0;
double ret = 0.; double ret = 0.;
if(A->Type != SCALAR){ if(A->Type != SCALAR) {
Message::Error("Non scalar argument for function 'AtIndex"); Message::Error("Non scalar argument for function 'AtIndex");
} }
else{ else {
index = (int)A->Val[0]; index = (int)A->Val[0];
if (index < 0 || index >= Fct->NbrParameters){ if(index < 0 || index >= Fct->NbrParameters) {
Message::Error("Wrong index in function 'AtIndex': %d not in [0,%d[", Message::Error("Wrong index in function 'AtIndex': %d not in [0,%d[",
index, Fct->NbrParameters); index, Fct->NbrParameters);
} }
else{ else {
ret = Fct->Para[index]; ret = Fct->Para[index];
} }
} }
V->Type = SCALAR; V->Type = SCALAR;
V->Val[0] = ret; V->Val[0] = ret;
if (Current.NbrHar != 1){ if(Current.NbrHar != 1) {
V->Val[MAX_DIM] = 0. ; V->Val[MAX_DIM] = 0.;
for (int k = 2 ; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar) ; k += 2) for(int k = 2; k < std::min(NBR_MAX_HARMONIC, Current.NbrHar); k += 2)
V->Val[MAX_DIM*k] = V->Val[MAX_DIM*(k+1)] = 0. ; V->Val[MAX_DIM * k] = V->Val[MAX_DIM * (k + 1)] = 0.;
} }
} }
 End of changes. 307 change blocks. 
887 lines changed or deleted 855 lines changed or added
To the C++ Programs section of the getdp source changes report or the top of this page
Mainly generated by pkgdiff using rfcdiff
Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)