#define RCSID "$Id: Get_Geometry.c,v 1.36 2006/02/26 00:42:54 geuzaine Exp $"
/*
* Copyright (C) 1997-2006 P. Dular, C. Geuzaine
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
* USA.
*
* Please report all bugs and problems to <getdp@geuz.org>.
*
* Contributor(s):
* Patrick Lefevre
*/
#if defined(HAVE_GSL)
#include <gsl/gsl_vector.h>
#include <gsl/gsl_linalg.h>
#else
#define NRANSI
#include "nrutil.h"
#endif
#include "GetDP.h"
#include "Get_Geometry.h"
#include "BF_Function.h"
#include "Numeric.h"
#include "Data_DefineE.h"
/* ------------------------------------------------------------------------ */
/* G e t _ N o d e s C o o r d i n a t e s O f E l e m e n t */
/* ------------------------------------------------------------------------ */
void Get_NodesCoordinatesOfElement(struct Element * Element) {
GetDP_Begin("Get_NodesCoordinatesOfElement");
if (Element->NumLastElementForNodesCoordinates != Element->Num) {
Element->NumLastElementForNodesCoordinates = Element->Num ;
Geo_GetNodesCoordinates
(Element->GeoElement->NbrNodes, Element->GeoElement->NumNodes,
Element->x, Element->y, Element->z) ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G e t _ G e o E l e m e n t */
/* ------------------------------------------------------------------------ */
void Get_BFGeoElement(struct Element * Element, double u, double v, double w) {
int i ;
GetDP_Begin("Get_BFGeoElement");
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
BF_Node (Element, i+1, u, v, w, &(Element->n[i])) ;
BF_GradNode(Element, i+1, u, v, w, Element->dndu[i]) ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G e t _ J a c o b i a n F u n c t i o n */
/* ------------------------------------------------------------------------ */
void * Get_JacobianFunction (int Type_Jacobian, int Type_Element,
int * Type_Dimension) {
GetDP_Begin("Get_JacobianFunction");
switch (Type_Jacobian) {
case JACOBIAN_VOL :
switch (Type_Element) {
case POINT :
*Type_Dimension = _0D ; GetDP_Return((void *)JacobianVol0D) ;
case LINE : case LINE_2 :
*Type_Dimension = _1D ; GetDP_Return((void *)JacobianVol1D) ;
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVol2D) ;
case TETRAHEDRON : case TETRAHEDRON_2 :
case HEXAHEDRON : case HEXAHEDRON_2 :
case PRISM : case PRISM_2 :
case PYRAMID : case PYRAMID_2 :
*Type_Dimension = _3D ; GetDP_Return((void *)JacobianVol3D) ;
default :
Msg(GERROR, "Unknown Jacobian Vol for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_SPH_SHELL :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolSphShell2D) ;
case TETRAHEDRON : case TETRAHEDRON_2 :
case HEXAHEDRON : case HEXAHEDRON_2 :
case PRISM : case PRISM_2 :
case PYRAMID : case PYRAMID_2 :
*Type_Dimension = _3D ; GetDP_Return((void *)JacobianVolSphShell3D) ;
default :
Msg(GERROR, "Unknown Jacobian VolSphShell for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_RECT_SHELL :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolRectShell2D) ;
case TETRAHEDRON : case TETRAHEDRON_2 :
case HEXAHEDRON : case HEXAHEDRON_2 :
case PRISM : case PRISM_2 :
case PYRAMID : case PYRAMID_2 :
*Type_Dimension = _3D ; GetDP_Return((void *)JacobianVolRectShell3D) ;
default :
Msg(GERROR, "Unknown Jacobian VolRectShell for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_PLPD_X :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolPlpdX2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolPlpdX for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxi2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxi for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI_SPH_SHELL :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSphShell2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxiSphShell for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI_RECT_SHELL :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiRectShell2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxiRectShell for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI_PLPD_X :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiPlpdX2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxiPlpdX for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI_SQU :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSqu2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxiSqu for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI_SQU_SPH_SHELL :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSquSphShell2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxiSquSphShell for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_VOL_AXI_SQU_RECT_SHELL :
switch (Type_Element) {
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVolAxiSquRectShell2D) ;
default :
Msg(GERROR, "Unknown Jacobian VolAxiSquRectShell for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_SUR :
switch (Type_Element) {
case POINT :
*Type_Dimension = _1D ; GetDP_Return((void *)JacobianVol0D) ;
case LINE : case LINE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianSur2D) ;
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
*Type_Dimension = _3D ; GetDP_Return((void *)JacobianSur3D) ;
default :
Msg(GERROR, "Unknown Jacobian Sur for element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_SUR_SPH_SHELL :
switch (Type_Element) {
case LINE : case LINE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianSurSphShell2D) ;
default :
Msg(GERROR, "Unknown Jacobian SurSphShell for element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_SUR_AXI :
switch (Type_Element) {
case LINE : case LINE_2 :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianSurAxi2D) ;
default :
Msg(GERROR, "Unknown Jacobian SurAxi for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
case JACOBIAN_LIN :
switch (Type_Element) {
case POINT :
*Type_Dimension = _2D ; GetDP_Return((void *)JacobianVol0D) ;
case LINE : case LINE_2 :
*Type_Dimension = _3D ; GetDP_Return((void *)JacobianLin3D) ;
default :
Msg(GERROR, "Unknown Jacobian Lin for Element Type (%s)",
Get_StringForDefine(Element_Type, Type_Element));
}
default :
Msg(GERROR, "Unknown Jacobian"); GetDP_Return(NULL) ;
}
}
/* ------------------------------------------------------------------------ */
/* G e o m e t r i c a l T r a n s f o r m a t i o n s */
/* ------------------------------------------------------------------------ */
double PlpdX2D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double CoorX, CoorY, A, B, R, theta, f ;
double DetJac ;
GetDP_Begin("PlpdX2D");
CoorX = CoorY = 0. ;
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
CoorX += Element->x[i] * Element->n[i] ;
CoorY += Element->y[i] * Element->n[i] ;
}
A = Element->JacobianCase->Para[0] ; B = Element->JacobianCase->Para[1] ;
R = CoorX ;
if ( (R > B+1.e-12*B) || (R < A-1.e-12*A) )
Msg(GERROR, "Bad parameters for unidirectional transformation Jacobian: "
"Rint=%g, Rext=%g, R=%g", A, B, R) ;
if (B == R) {
Jac->c11 = 1. ; Jac->c12 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 1. ;
GetDP_Return(1.) ;
}
f = (A*(B-A)) / (R*(B-R)) ;
theta = (B-2.*R) / (B-R) ;
Jac->c11 = f * (1.- theta) ; Jac->c12 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 1. ;
DetJac = f*( 1.-theta) ;
GetDP_Return(DetJac) ;
}
double power(double x, double y) {
if (y == 1.0) return (x);
else if (y == 2.0) return (x*x);
else if (y == 0.5) return (sqrt(x));
else return (pow(x,y));
}
double Transformation (int Dim, int Type, struct Element * Element, MATRIX3x3 * Jac) {
int i, Axis = 0 ;
double X = 0., Y = 0., Z = 0. ;
double p = 1., L= 0. ;
double Cx = 0., Cy = 0., Cz = 0., A = 0., B = 0., R = 0. ;
double theta, XR, YR, ZR, f, dRdx = 0., dRdy = 0., dRdz = 0. ;
double DetJac ;
/*
A = interior radius
B = exterior radius
Cx, Cy, Cz = coord of centre
Axis = direction of the transformation
p = exponant
1/L = f(B)
*/
GetDP_Begin("Transformation");
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++) {
X += Element->x[i] * Element->n[i] ;
Y += Element->y[i] * Element->n[i] ;
Z += Element->z[i] * Element->n[i] ;
}
if(Element->JacobianCase->NbrParameters >= 2){
A = Element->JacobianCase->Para[0] ;
B = Element->JacobianCase->Para[1] ;
}
else
Msg(GERROR, "Missing interior and/or exterior radius for transformation Jacobian");
if(Type == JACOBIAN_RECT){
if(Element->JacobianCase->NbrParameters >= 3)
Axis = (int)Element->JacobianCase->Para[2] ;
if(Element->JacobianCase->NbrParameters >= 6){
Cx = Element->JacobianCase->Para[3] ;
Cy = Element->JacobianCase->Para[4] ;
Cz = Element->JacobianCase->Para[5] ;
}
if(Element->JacobianCase->NbrParameters >= 7)
p = Element->JacobianCase->Para[6] ;
if(Element->JacobianCase->NbrParameters >= 8)
L = Element->JacobianCase->Para[7] ;
if(Element->JacobianCase->NbrParameters >= 9){
Msg(DIRECT, ERROR_STR "Too many parameters for rectangular transformation Jacobian");
Msg(GERROR, "Valid parameters: Radius1 Radius2 Axis CenterX CenterY CenterZ Power 1/Infinity");
}
}
else if(Type == JACOBIAN_SPH){
if(Element->JacobianCase->NbrParameters >= 5){
Cx = Element->JacobianCase->Para[2] ;
Cy = Element->JacobianCase->Para[3] ;
Cz = Element->JacobianCase->Para[4] ;
}
if(Element->JacobianCase->NbrParameters >= 6)
p = Element->JacobianCase->Para[5] ;
if(Element->JacobianCase->NbrParameters >= 7)
L = Element->JacobianCase->Para[6] ;
if(Element->JacobianCase->NbrParameters >= 8){
Msg(DIRECT, ERROR_STR "Too many parameters for spherical transformation Jacobian");
Msg(GERROR, "Valid parameters: Radius1 Radius2 CenterX CenterY CenterZ Power 1/Infinity");
}
}
else
Msg(GERROR, "Unknown type of transformation Jacobian");
if(L) B = (B*B-A*A*L)/(B-A*L);
if(Type == JACOBIAN_SPH){
R = sqrt( DSQU(X-Cx) + DSQU(Y-Cy) + DSQU(Z-Cz) ) ;
dRdx = (X-Cx)/R ;
dRdy = (Y-Cy)/R ;
dRdz = (Z-Cz)/R ;
}
else{
switch(Axis) {
case 1: R = (X-Cx) ; dRdx =1.0 ; dRdy =0.0 ; dRdz =0.0 ; break;
case 2: R = (Y-Cy) ; dRdx =0.0 ; dRdy =1.0 ; dRdz =0.0 ; break;
case 3: R = (Z-Cz) ; dRdx =0.0 ; dRdy =0.0 ; dRdz =1.0 ; break;
default: Msg(GERROR, "Bad axis specification: 1 for X, 2 for Y and 3 for Z");
}
}
if ( (fabs(R) > fabs(B) + 1.e-12*fabs(B)) ||
(fabs(R) < fabs(A) - 1.e-2*fabs(A)) )
Msg(GERROR, "Bad parameters for transformation Jacobian: %g not in [%g,%g]", R, A, B) ;
if (B == R) {
Jac->c11 = 1. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 1. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
GetDP_Return(1.) ;
}
f = power((A*(B-A))/(R*(B-R)), p) ;
theta = p * (B-2.*R) / (B-R) ;
XR = (X-Cx)/R ;
YR = (Y-Cy)/R ;
ZR = (Z-Cz)/R ;
Jac->c11 = f * (1.0 - theta * XR * dRdx) ;
Jac->c12 = f * ( - theta * XR * dRdy) ;
Jac->c13 = f * ( - theta * XR * dRdz) ;
Jac->c21 = f * ( - theta * YR * dRdx) ;
Jac->c22 = f * (1.0 - theta * YR * dRdy) ;
Jac->c23 = f * ( - theta * YR * dRdz) ;
Jac->c31 = f * ( - theta * ZR * dRdx) ;
Jac->c32 = f * ( - theta * ZR * dRdy) ;
Jac->c33 = f * (1.0 - theta * ZR * dRdz) ;
DetJac = f * f * (1.0 - theta) ;
GetDP_Return(DetJac) ;
}
/* ------------------------------------------------------------------------ */
/* J a c o b i a n V o l */
/* ------------------------------------------------------------------------ */
/* 0D */
double JacobianVol0D (struct Element * Element, MATRIX3x3 * Jac) {
GetDP_Begin("JacobianVol0D");
Jac->c11 = 1. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 1. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
GetDP_Return(1.) ;
}
/* 1D */
double JacobianVol1D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianVol1D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 1. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
}
DetJac = Jac->c11 ;
GetDP_Return(DetJac) ;
}
/* 2D */
double JacobianVol2D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianVol2D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 0. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
Jac->c21 += Element->x[i] * Element->dndu[i][1] ;
Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
Jac->c22 += Element->y[i] * Element->dndu[i][1] ;
}
DetJac = Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21 ;
GetDP_Return(DetJac) ;
}
double JacobianVolSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolSphShell2D");
DetJac1 = JacobianVol2D (Element, &Jac1) ;
DetJac2 = Transformation(_2D, JACOBIAN_SPH, Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianVolRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolRectShell2D");
DetJac1 = JacobianVol2D (Element, &Jac1) ;
DetJac2 = Transformation (_2D, JACOBIAN_RECT, Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianVolPlpdX2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolPlpdX2D");
DetJac1 = JacobianVol2D (Element, &Jac1) ;
DetJac2 = PlpdX2D (Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
GetDP_Return(DetJac1 * DetJac2) ;
}
/* 2D Axi (Attention, l'axe doit etre x=z=0) */
double JacobianVolAxi2D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double s = 0., DetJac ;
GetDP_Begin("JacobianVolAxi2D");
DetJac = JacobianVol2D(Element, Jac) ;
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
s += Element->x[i] * Element->n[i] ;
/* Warning! For evaluations on the symmetry axis */
if (s==0.0) {
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
s += Element->x[i] ;
s /= (double)Element->GeoElement->NbrNodes ;
}
Jac->c33 = s ;
GetDP_Return(DetJac * Jac->c33) ;
}
double JacobianVolAxiSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolAxiSphShell2D");
DetJac1 = JacobianVolAxi2D (Element, &Jac1) ;
DetJac2 = Transformation (_2D, JACOBIAN_SPH, Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = Jac1.c33 ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianVolAxiRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolAxiRectShell2D");
DetJac1 = JacobianVolAxi2D (Element, &Jac1) ;
DetJac2 = Transformation (_2D, JACOBIAN_RECT, Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = Jac1.c33 ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianVolAxiPlpdX2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolAxiPlpdX2D");
DetJac1 = JacobianVolAxi2D (Element, &Jac1) ;
DetJac2 = PlpdX2D (Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = Jac1.c33 ;
GetDP_Return(DetJac1 * DetJac2) ;
}
/* 2D Axi avec transformation quadratique (Attention, l'axe doit etre x=z=0) */
double JacobianVolAxiSqu2D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double s = 0., r, DetJac ;
GetDP_Begin("JacobianVolAxiSqu2D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 0. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
s += DSQU(Element->x[i]) * Element->n[i] ;
/* Warning! For evaluations on the symmetry axis */
if (s==0.0) {
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
s += Element->x[i] * Element->x[i] ;
s /= (double)Element->GeoElement->NbrNodes ;
}
r = sqrt(s);
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += 0.5/r * DSQU(Element->x[i]) * Element->dndu[i][0] ;
Jac->c21 += 0.5/r * DSQU(Element->x[i]) * Element->dndu[i][1] ;
Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
Jac->c22 += Element->y[i] * Element->dndu[i][1] ;
}
Jac->c33 = r ;
DetJac = (Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21) * Jac->c33 ;
GetDP_Return(DetJac) ;
}
double JacobianVolAxiSquSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolAxiSquSphShell2D");
DetJac1 = JacobianVolAxiSqu2D(Element, &Jac1) ;
DetJac2 = Transformation (_2D, JACOBIAN_SPH, Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = Jac1.c33 ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianVolAxiSquRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolAxiSquRectShell2D");
DetJac1 = JacobianVolAxiSqu2D(Element, &Jac1) ;
DetJac2 = Transformation (_2D, JACOBIAN_RECT, Element, &Jac2) ;
Get_ProductMatrix( _2D, &Jac1, &Jac2, Jac) ;
Jac->c13 = 0. ;
Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = Jac1.c33 ;
GetDP_Return(DetJac1 * DetJac2) ;
}
/* 3D */
double JacobianVol3D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianVol3D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 0. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 0. ;
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
Jac->c21 += Element->x[i] * Element->dndu[i][1] ;
Jac->c31 += Element->x[i] * Element->dndu[i][2] ;
Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
Jac->c22 += Element->y[i] * Element->dndu[i][1] ;
Jac->c32 += Element->y[i] * Element->dndu[i][2] ;
Jac->c13 += Element->z[i] * Element->dndu[i][0] ;
Jac->c23 += Element->z[i] * Element->dndu[i][1] ;
Jac->c33 += Element->z[i] * Element->dndu[i][2] ;
}
DetJac = Jac->c11 * Jac->c22 * Jac->c33 + Jac->c13 * Jac->c21 * Jac->c32
+ Jac->c12 * Jac->c23 * Jac->c31 - Jac->c13 * Jac->c22 * Jac->c31
- Jac->c11 * Jac->c23 * Jac->c32 - Jac->c12 * Jac->c21 * Jac->c33 ;
GetDP_Return(DetJac) ;
}
double JacobianVolSphShell3D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolShell3D");
DetJac1 = JacobianVol3D (Element, &Jac1) ;
DetJac2 = Transformation(_3D, JACOBIAN_SPH, Element, &Jac2) ;
Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianVolRectShell3D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianVolRectShell3D");
DetJac1 = JacobianVol3D (Element, &Jac1) ;
DetJac2 = Transformation (_3D, JACOBIAN_RECT, Element, &Jac2) ;
Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;
GetDP_Return(DetJac1 * DetJac2) ;
}
/* ------------------------------------------------------------------------ */
/* J a c o b i a n S u r */
/* ------------------------------------------------------------------------ */
double JacobianSur2D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianSur2D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 0. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 1. ;
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
}
DetJac = HYPOT(Jac->c11, Jac->c12) ;
GetDP_Return(DetJac) ;
}
double JacobianSurSphShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianSurSphShell2D");
DetJac1 = JacobianSur2D (Element, &Jac1) ;
DetJac2 = Transformation(_2D, JACOBIAN_SPH, Element, &Jac2) ;
Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianSurRectShell2D (struct Element * Element, MATRIX3x3 * Jac) {
MATRIX3x3 Jac1, Jac2 ;
double DetJac1, DetJac2 ;
GetDP_Begin("JacobianSurRectShell2D");
DetJac1 = JacobianSur2D (Element, &Jac1) ;
DetJac2 = Transformation(_2D, JACOBIAN_RECT, Element, &Jac2) ;
Get_ProductMatrix( _3D, &Jac1, &Jac2, Jac) ;
GetDP_Return(DetJac1 * DetJac2) ;
}
double JacobianSurAxi2D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianSurAxi2D");
DetJac = JacobianSur2D(Element, Jac) ;
Jac->c33 = 0. ;
for (i = 0 ; i < Element->GeoElement->NbrNodes ; i++)
Jac->c33 += Element->x[i] * Element->n[i] ;
GetDP_Return(DetJac * Jac->c33) ;
}
double JacobianSur3D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianSur3D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 0. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 0. ;
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
Jac->c21 += Element->x[i] * Element->dndu[i][1] ;
Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
Jac->c22 += Element->y[i] * Element->dndu[i][1] ;
Jac->c13 += Element->z[i] * Element->dndu[i][0] ;
Jac->c23 += Element->z[i] * Element->dndu[i][1] ;
}
DetJac = sqrt( DSQU(Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21)
+ DSQU(Jac->c13 * Jac->c21 - Jac->c11 * Jac->c23)
+ DSQU(Jac->c12 * Jac->c23 - Jac->c13 * Jac->c22) ) ;
GetDP_Return(DetJac) ;
}
/* ------------------------------------------------------------------------ */
/* J a c o b i a n L i n */
/* ------------------------------------------------------------------------ */
double JacobianLin3D (struct Element * Element, MATRIX3x3 * Jac) {
int i ;
double DetJac ;
GetDP_Begin("JacobianLin3D");
Jac->c11 = 0. ; Jac->c12 = 0. ; Jac->c13 = 0. ;
Jac->c21 = 0. ; Jac->c22 = 0. ; Jac->c23 = 0. ;
Jac->c31 = 0. ; Jac->c32 = 0. ; Jac->c33 = 0. ;
for ( i = 0 ; i < Element->GeoElement->NbrNodes ; i++ ) {
Jac->c11 += Element->x[i] * Element->dndu[i][0] ;
Jac->c12 += Element->y[i] * Element->dndu[i][0] ;
Jac->c13 += Element->z[i] * Element->dndu[i][0] ;
}
DetJac = sqrt(DSQU(Jac->c11)+DSQU(Jac->c12)+DSQU(Jac->c13)) ;
GetDP_Return(DetJac) ;
}
/* ------------------------------------------------------------------------ */
/* G e t _ I n v e r s e M a t r i x */
/* ------------------------------------------------------------------------ */
void Get_InverseSingularMatrix(MATRIX3x3 * Mat, MATRIX3x3 * InvMat) {
double T[3][3] = {{0.,0.,0.},{0.,0.,0.},{0.,0.,0.}};
int i, j, k;
#if defined(HAVE_GSL)
gsl_matrix *M, *V;
gsl_vector *W, *TMPVEC;
double ww;
#else
double **M, **V, *W;
#endif
GetDP_Begin("Get_InverseSingularMatrix");
InvMat->c11 = InvMat->c12 = InvMat->c13 = 0.0;
InvMat->c21 = InvMat->c22 = InvMat->c23 = 0.0;
InvMat->c31 = InvMat->c32 = InvMat->c33 = 0.0;
#if defined(HAVE_GSL)
M = gsl_matrix_alloc(3, 3);
V = gsl_matrix_alloc(3, 3);
W = gsl_vector_alloc(3);
TMPVEC = gsl_vector_alloc(3);
gsl_matrix_set(M, 0, 0, Mat->c11);
gsl_matrix_set(M, 0, 1, Mat->c12);
gsl_matrix_set(M, 0, 2, Mat->c13);
gsl_matrix_set(M, 1, 0, Mat->c21);
gsl_matrix_set(M, 1, 1, Mat->c22);
gsl_matrix_set(M, 1, 2, Mat->c23);
gsl_matrix_set(M, 2, 0, Mat->c31);
gsl_matrix_set(M, 2, 1, Mat->c32);
gsl_matrix_set(M, 2, 2, Mat->c33);
gsl_linalg_SV_decomp(M, V, W, TMPVEC);
for(i = 0; i < 3; i++) {
for(j = 0; j < 3; j++) {
ww = gsl_vector_get(W, i);
if(fabs(ww) > 1.e-16) { /* singular value precision */
T[i][j] += gsl_matrix_get(M, j, i) / ww;
}
}
}
for(k=0 ; k<3 ; k++){
InvMat->c11 += gsl_matrix_get(V, 0, k) * T[k][0] ;
InvMat->c12 += gsl_matrix_get(V, 0, k) * T[k][1] ;
InvMat->c13 += gsl_matrix_get(V, 0, k) * T[k][2] ;
InvMat->c21 += gsl_matrix_get(V, 1, k) * T[k][0] ;
InvMat->c22 += gsl_matrix_get(V, 1, k) * T[k][1] ;
InvMat->c23 += gsl_matrix_get(V, 1, k) * T[k][2] ;
InvMat->c31 += gsl_matrix_get(V, 2, k) * T[k][0] ;
InvMat->c32 += gsl_matrix_get(V, 2, k) * T[k][1] ;
InvMat->c33 += gsl_matrix_get(V, 2, k) * T[k][2] ;
}
gsl_matrix_free(M);
gsl_matrix_free(V);
gsl_vector_free(W);
gsl_vector_free(TMPVEC);
#else
M = dmatrix(1,3,1,3);
V = dmatrix(1,3,1,3);
W = dvector(1,3);
M[1][1] = Mat->c11 ; M[1][2] = Mat->c12 ; M[1][3] = Mat->c13 ;
M[2][1] = Mat->c21 ; M[2][2] = Mat->c22 ; M[2][3] = Mat->c23 ;
M[3][1] = Mat->c31 ; M[3][2] = Mat->c32 ; M[3][3] = Mat->c33 ;
dsvdcmp(M, 3, 3, W, V);
/* cf. Numerical Recipes in C, p. 62 */
for(i=1 ; i<=3 ; i++)
for(j=1 ; j<=3 ; j++)
if(fabs(W[i]) > 1.e-16) /* precision */
T[i-1][j-1] += M[j][i] / W[i] ;
for(k=1 ; k<=3 ; k++){
InvMat->c11 += V[1][k] * T[k-1][0] ;
InvMat->c12 += V[1][k] * T[k-1][1] ;
InvMat->c13 += V[1][k] * T[k-1][2] ;
InvMat->c21 += V[2][k] * T[k-1][0] ;
InvMat->c22 += V[2][k] * T[k-1][1] ;
InvMat->c23 += V[2][k] * T[k-1][2] ;
InvMat->c31 += V[3][k] * T[k-1][0] ;
InvMat->c32 += V[3][k] * T[k-1][1] ;
InvMat->c33 += V[3][k] * T[k-1][2] ;
}
free_dmatrix(M,1,3,1,3);
free_dmatrix(V,1,3,1,3);
free_dvector(W,1,3);
#endif
GetDP_End ;
}
void Get_InverseMatrix(int Type_Dimension, int Type_Element, double DetMat,
MATRIX3x3 * Mat, MATRIX3x3 * InvMat) {
GetDP_Begin("Get_InverseMatrix");
switch (Type_Dimension) {
case _0D :
InvMat->c11 = InvMat->c22 = InvMat->c33 = 1. ;
InvMat->c12 = InvMat->c21 = 0. ;
InvMat->c13 = InvMat->c31 = 0. ;
InvMat->c23 = InvMat->c32 = 0. ;
break ;
case _1D :
InvMat->c11 = 1. / Mat->c11 ;
InvMat->c22 = 1. / Mat->c22 ;
InvMat->c33 = 1. / Mat->c33 ;
InvMat->c12 = InvMat->c21 = 0. ;
InvMat->c13 = InvMat->c31 = 0. ;
InvMat->c23 = InvMat->c32 = 0. ;
break ;
case _2D :
switch(Type_Element){
case TRIANGLE : case TRIANGLE_2 :
case QUADRANGLE : case QUADRANGLE_2 :
if(!DetMat) Msg(GERROR, "Null determinant in 'Get_InverseMatrix'");
InvMat->c11 = Mat->c22 * Mat->c33 / DetMat ;
InvMat->c21 = - Mat->c21 * Mat->c33 / DetMat ;
InvMat->c12 = - Mat->c12 * Mat->c33 / DetMat ;
InvMat->c22 = Mat->c11 * Mat->c33 / DetMat ;
InvMat->c13 = InvMat->c23 = InvMat->c31 = InvMat->c32 = 0. ;
InvMat->c33 = 1. / Mat->c33 ;
break ;
default :
Get_InverseSingularMatrix(Mat, InvMat);
break;
}
break;
case _3D :
switch(Type_Element){
case TETRAHEDRON : case TETRAHEDRON_2 :
case HEXAHEDRON : case HEXAHEDRON_2 :
case PRISM : case PRISM_2 :
case PYRAMID : case PYRAMID_2 :
if(!DetMat) Msg(GERROR, "Null determinant in 'Get_InverseMatrix'");
InvMat->c11 = ( Mat->c22 * Mat->c33 - Mat->c23 * Mat->c32 ) / DetMat ;
InvMat->c21 = -( Mat->c21 * Mat->c33 - Mat->c23 * Mat->c31 ) / DetMat ;
InvMat->c31 = ( Mat->c21 * Mat->c32 - Mat->c22 * Mat->c31 ) / DetMat ;
InvMat->c12 = -( Mat->c12 * Mat->c33 - Mat->c13 * Mat->c32 ) / DetMat ;
InvMat->c22 = ( Mat->c11 * Mat->c33 - Mat->c13 * Mat->c31 ) / DetMat ;
InvMat->c32 = -( Mat->c11 * Mat->c32 - Mat->c12 * Mat->c31 ) / DetMat ;
InvMat->c13 = ( Mat->c12 * Mat->c23 - Mat->c13 * Mat->c22 ) / DetMat ;
InvMat->c23 = -( Mat->c11 * Mat->c23 - Mat->c13 * Mat->c21 ) / DetMat ;
InvMat->c33 = ( Mat->c11 * Mat->c22 - Mat->c12 * Mat->c21 ) / DetMat ;
break;
default :
Get_InverseSingularMatrix(Mat, InvMat);
break;
}
break;
default :
Msg(GERROR, "Wrong dimension in 'Get_InverseMatrix'");
break ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G e t _ P r o d u c t M a t r i x */
/* ------------------------------------------------------------------------ */
void Get_ProductMatrix(int Type_Dimension,
MATRIX3x3 * A, MATRIX3x3 * B, MATRIX3x3 * AB) {
GetDP_Begin("Get_ProductMatrix");
switch (Type_Dimension) {
case _2D :
AB->c11 = A->c11 * B->c11 + A->c12 * B->c21 ;
AB->c12 = A->c11 * B->c12 + A->c12 * B->c22 ;
AB->c21 = A->c21 * B->c11 + A->c22 * B->c21 ;
AB->c22 = A->c21 * B->c12 + A->c22 * B->c22 ;
break ;
case _3D :
AB->c11 = A->c11 * B->c11 + A->c12 * B->c21 + A->c13 * B->c31 ;
AB->c12 = A->c11 * B->c12 + A->c12 * B->c22 + A->c13 * B->c32 ;
AB->c13 = A->c11 * B->c13 + A->c12 * B->c23 + A->c13 * B->c33 ;
AB->c21 = A->c21 * B->c11 + A->c22 * B->c21 + A->c23 * B->c31 ;
AB->c22 = A->c21 * B->c12 + A->c22 * B->c22 + A->c23 * B->c32 ;
AB->c23 = A->c21 * B->c13 + A->c22 * B->c23 + A->c23 * B->c33 ;
AB->c31 = A->c31 * B->c11 + A->c32 * B->c21 + A->c33 * B->c31 ;
AB->c32 = A->c31 * B->c12 + A->c32 * B->c22 + A->c33 * B->c32 ;
AB->c33 = A->c31 * B->c13 + A->c32 * B->c23 + A->c33 * B->c33 ;
break ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* G e t _ C h a n g e O f C o o r d i n a t e s */
/* ------------------------------------------------------------------------ */
void *Get_ChangeOfCoordinates(int Flag_ChangeCoord, int Type_Form) {
GetDP_Begin("Get_ChangeOfCoordinates");
switch (Type_Form) {
case SCALAR :
case FORM0 :
GetDP_Return((void *)ChangeOfCoord_No1) ;
case FORM1 :
GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form1 :
(void *)ChangeOfCoord_No123) ;
case FORM2 :
GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form2 :
(void *)ChangeOfCoord_No123) ;
case FORM3 : case FORM3P :
GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form3 :
(void *)ChangeOfCoord_No1) ;
case FORM1P :
GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form1P :
(void *)ChangeOfCoord_No123) ;
case FORM2P :
GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form2P :
(void *)ChangeOfCoord_No123) ;
case VECTOR :
GetDP_Return((void *)ChangeOfCoord_No123) ;
case FORM1S :
GetDP_Return((Flag_ChangeCoord) ? (void *)ChangeOfCoord_Form1S :
(void *)ChangeOfCoord_No123) ;
default :
Msg(GERROR, "Unknown type of field (%s)",
Get_StringForDefine(Field_Type, Type_Form)) ;
GetDP_Return(NULL) ;
}
}
/* ------------------------------------------------------------------------ */
/* C h a n g e O f C o o r d _ X X X */
/* ------------------------------------------------------------------------ */
void ChangeOfCoord_No1(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_No1");
vBFx[0] = vBFu[0] ;
GetDP_End ;
}
void ChangeOfCoord_No123(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_No123");
vBFx[0] = vBFu[0] ; vBFx[1] = vBFu[1] ; vBFx[2] = vBFu[2] ;
GetDP_End ;
}
void ChangeOfCoord_Form1(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_Form1");
vBFx[0] = vBFu[0] * Element->InvJac.c11 + vBFu[1] * Element->InvJac.c12
+ vBFu[2] * Element->InvJac.c13 ;
vBFx[1] = vBFu[0] * Element->InvJac.c21 + vBFu[1] * Element->InvJac.c22
+ vBFu[2] * Element->InvJac.c23 ;
vBFx[2] = vBFu[0] * Element->InvJac.c31 + vBFu[1] * Element->InvJac.c32
+ vBFu[2] * Element->InvJac.c33 ;
GetDP_End ;
}
void ChangeOfCoord_Form2(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_Form2");
if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form2'");
vBFx[0] = (vBFu[0] * Element->Jac.c11 + vBFu[1] * Element->Jac.c21
+ vBFu[2] * Element->Jac.c31) / Element->DetJac ;
vBFx[1] = (vBFu[0] * Element->Jac.c12 + vBFu[1] * Element->Jac.c22
+ vBFu[2] * Element->Jac.c32) / Element->DetJac ;
vBFx[2] = (vBFu[0] * Element->Jac.c13 + vBFu[1] * Element->Jac.c23
+ vBFu[2] * Element->Jac.c33) / Element->DetJac ;
GetDP_End ;
}
void ChangeOfCoord_Form3(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_Form3");
if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form3'");
vBFx[0] = vBFu[0] / Element->DetJac ;
GetDP_End ;
}
/* Form1P, 2P, 1S : valid in 2D only ! */
void ChangeOfCoord_Form1P(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_Form1P");
vBFx[0] = 0. ;
vBFx[1] = 0. ;
vBFx[2] = vBFu[2] / Element->Jac.c33 ; /* ... * Element->InvJac.c33 */
GetDP_End ;
}
void ChangeOfCoord_Form2P(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_Form2P");
if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form2P'");
vBFx[0] = (vBFu[0] * Element->Jac.c11 + vBFu[1] * Element->Jac.c21)
/ Element->DetJac ;
vBFx[1] = (vBFu[0] * Element->Jac.c12 + vBFu[1] * Element->Jac.c22)
/ Element->DetJac ;
vBFx[2] = 0. ;
GetDP_End ;
}
void ChangeOfCoord_Form1S(struct Element * Element,
double vBFu[], double vBFx[]) {
GetDP_Begin("ChangeOfCoord_Form1S");
if(!Element->DetJac) Msg(GERROR, "Null determinant in 'ChangeOfCoord_Form1S'");
vBFx[0] = 0. ;
vBFx[1] = 0. ;
vBFx[2] = vBFu[0] / Element->DetJac ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* C a l _ P r o d u c t X X X */
/* ------------------------------------------------------------------------ */
double Cal_Product123(double v1[], double v2[]) {
return v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2] ;
}
double Cal_Product12 (double v1[], double v2[]) {
return v1[0]*v2[0] + v1[1]*v2[1] ;
}
double Cal_Product3 (double v1[], double v2[]) {
return v1[2]*v2[2] ;
}
double Cal_Product1 (double v1[], double v2[]) {
return v1[0]*v2[0] ;
}
syntax highlighted by Code2HTML, v. 0.9.1