#define RCSID "$Id: F_Misc.c,v 1.32 2006/03/02 22:04:12 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):
* Tuan Ledinh
*/
#include "GetDP.h"
#include "Data_DefineE.h"
#include "F_Function.h"
#include "GeoData.h"
#include "Get_Geometry.h"
#include "Cal_Value.h"
#include "CurrentData.h"
#include "Numeric.h"
#include "Tools.h"
#if !defined(HAVE_GSL)
#define NRANSI
#include "nrutil.h" /* pour Tuan */
#endif
/* ------------------------------------------------------------------------ */
/* Warning: the pointers A and V can be identical. You must */
/* use temporary variables in your computations: you can only */
/* affect to V at the very last time (when you're sure you will */
/* not use A any more). */
/* ------------------------------------------------------------------------ */
#define F_ARG struct Function * Fct, struct Value * A, struct Value * V
/* ------------------------------------------------------------------------ */
/* Printf */
/* ------------------------------------------------------------------------ */
void F_Printf (F_ARG) {
GetDP_Begin("F_Printf");
V = A ; /* Attention: ne sert a rien!?! */
Print_Value(A) ;
printf("\n") ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* Computes the normal to an element */
/* ------------------------------------------------------------------------ */
void F_Normal(F_ARG) {
int k ;
GetDP_Begin("F_Normal");
if(!Current.Element || Current.Element->Num == NO_ELEMENT)
Msg(GERROR, "No element on which to compute 'F_Normal'");
Geo_CreateNormal(Current.Element->Type,
Current.Element->x,
Current.Element->y,
Current.Element->z,
V->Val);
if (Current.NbrHar != 1) {
V->Val[MAX_DIM] = 0. ;
V->Val[MAX_DIM+1] = 0. ;
V->Val[MAX_DIM+2] = 0. ;
for (k = 2 ; k < Current.NbrHar ; k += 2) {
V->Val[MAX_DIM* k ] = V->Val[0] ;
V->Val[MAX_DIM* k +1] = V->Val[1] ;
V->Val[MAX_DIM* k +2] = V->Val[2] ;
V->Val[MAX_DIM*(k+1) ] = 0. ;
V->Val[MAX_DIM*(k+1)+1] = 0. ;
V->Val[MAX_DIM*(k+1)+2] = 0. ;
}
}
V->Type = VECTOR ;
GetDP_End ;
}
void F_NormalSource(F_ARG) {
int k ;
GetDP_Begin("F_NormalSource");
if(!Current.ElementSource || Current.ElementSource->Num == NO_ELEMENT)
Msg(GERROR, "No element on which to compute 'F_NormalSource'");
Geo_CreateNormal(Current.ElementSource->Type,
Current.ElementSource->x,
Current.ElementSource->y,
Current.ElementSource->z,
V->Val);
if (Current.NbrHar != 1) {
V->Val[MAX_DIM] = 0. ;
V->Val[MAX_DIM+1] = 0. ;
V->Val[MAX_DIM+2] = 0. ;
for (k = 2 ; k < Current.NbrHar ; k += 2) {
V->Val[MAX_DIM* k ] = V->Val[0] ;
V->Val[MAX_DIM* k +1] = V->Val[1] ;
V->Val[MAX_DIM* k +2] = V->Val[2] ;
V->Val[MAX_DIM*(k+1) ] = 0. ;
V->Val[MAX_DIM*(k+1)+1] = 0. ;
V->Val[MAX_DIM*(k+1)+2] = 0. ;
}
}
V->Type = VECTOR ;
GetDP_End ;
}
void F_Tangent(F_ARG) {
int k ;
double tx, ty, tz, norm ;
GetDP_Begin("F_Tangent");
if(!Current.Element || Current.Element->Num == NO_ELEMENT)
Msg(GERROR, "No element on which to compute 'F_Tangent'");
switch (Current.Element->Type) {
case LINE :
tx = Current.Element->x[1] - Current.Element->x[0] ;
ty = Current.Element->y[1] - Current.Element->y[0] ;
tz = Current.Element->z[1] - Current.Element->z[0] ;
norm = sqrt(DSQU(tx)+DSQU(ty)+DSQU(tz)) ;
V->Val[0] = tx/norm ;
V->Val[1] = ty/norm ;
V->Val[2] = tz/norm ;
break ;
default :
Msg(GERROR, "Function 'Tangent' only valid for Line Elements");
}
if (Current.NbrHar != 1) {
V->Val[MAX_DIM] = 0. ;
V->Val[MAX_DIM+1] = 0. ;
V->Val[MAX_DIM+2] = 0. ;
for (k = 2 ; k < Current.NbrHar ; k += 2) {
V->Val[MAX_DIM* k ] = V->Val[0] ;
V->Val[MAX_DIM* k +1] = V->Val[1] ;
V->Val[MAX_DIM* k +2] = V->Val[2] ;
V->Val[MAX_DIM*(k+1) ] = 0. ;
V->Val[MAX_DIM*(k+1)+1] = 0. ;
V->Val[MAX_DIM*(k+1)+2] = 0. ;
}
}
V->Type = VECTOR ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* Comparison */
/* ------------------------------------------------------------------------ */
void F_CompElementNum (F_ARG) {
GetDP_Begin("F_CompElementNum");
if(!Current.Element || !Current.ElementSource)
Msg(GERROR, "Uninitialized Element in 'F_CompElementNum'");
V->Type = SCALAR ;
V->Val[0] = (Current.Element->Num == Current.ElementSource->Num) ;
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* Volume */
/* ------------------------------------------------------------------------ */
void F_ElementVol (F_ARG) {
int k;
double Vol = 0.;
MATRIX3x3 Jac;
/* It would be more efficient to compute the volumes directly from
the node coordinates, but I'm lazy... */
Get_NodesCoordinatesOfElement(Current.Element) ;
Get_BFGeoElement(Current.Element, Current.u, Current.v, Current.w) ;
/* we use the most general case (3D embedding) */
switch(Current.Element->Type){
case LINE:
Vol = 2. * JacobianLin3D(Current.Element,&Jac);
break;
case TRIANGLE:
Vol = 0.5 * JacobianSur3D(Current.Element,&Jac) ;
break;
case QUADRANGLE:
Vol = 4. * JacobianSur3D(Current.Element,&Jac) ;
break;
case TETRAHEDRON:
Vol = 1./6. * JacobianVol3D(Current.Element,&Jac) ;
break;
case HEXAHEDRON:
Vol = 8. * JacobianVol3D(Current.Element,&Jac) ;
break;
case PRISM:
Vol = JacobianVol3D(Current.Element,&Jac) ;
break;
default :
Msg(GERROR, "F_ElementVol not implemented for %s",
Get_StringForDefine(Element_Type, Current.Element->Type));
}
V->Type = SCALAR ;
V->Val[0] = fabs(Vol);
for (k = 2 ; k < Current.NbrHar ; k += 2)
V->Val[MAX_DIM* k] = V->Val[0] ;
}
/* ------------------------------------------------------------------------ */
/* SurfaceArea */
/* ------------------------------------------------------------------------ */
void F_SurfaceArea (F_ARG) {
struct Element Element ;
List_T * InitialList_L;
int Index_Region, Nbr_Element, i_Element ;
double Val_Surface ;
double c11, c21, c12, c22, DetJac ;
int i, k ;
GetDP_Begin("F_SurfaceArea");
if (!Fct->Active) {
Fct->Active = (struct FunctionActive *)Malloc(sizeof(struct FunctionActive)) ;
if (Fct->NbrParameters == 1) {
Index_Region = (int)(Fct->Para[0]) ;
InitialList_L = List_Create(1,1,sizeof(int));
List_Reset(InitialList_L);
List_Add(InitialList_L,&Index_Region);
/*
InitialList_L = ((struct Group *)
List_Pointer(Problem_S.Group, Index_Region))->InitialList ;
*/
}
else {
Index_Region = -1 ;
InitialList_L = NULL ;
}
Val_Surface = 0. ;
Nbr_Element = Geo_GetNbrGeoElements() ;
for (i_Element = 0 ; i_Element < Nbr_Element; i_Element++) {
Element.GeoElement = Geo_GetGeoElement(i_Element) ;
if ((InitialList_L &&
List_Search(InitialList_L, &(Element.GeoElement->Region), fcmp_int)) ||
(!InitialList_L && Element.GeoElement->Region == Current.Region)) {
Element.Num = Element.GeoElement->Num ;
Element.Type = Element.GeoElement->Type ;
if (Element.Type == TRIANGLE || Element.Type == QUADRANGLE) {
Get_NodesCoordinatesOfElement(&Element) ;
Get_BFGeoElement(&Element, 0., 0., 0.) ;
c11 = c21 = c12 = c22 = 0. ;
for ( i = 0 ; i < Element.GeoElement->NbrNodes ; i++ ) {
c11 += Element.x[i] * Element.dndu[i][0] ;
c21 += Element.x[i] * Element.dndu[i][1] ;
c12 += Element.y[i] * Element.dndu[i][0] ;
c22 += Element.y[i] * Element.dndu[i][1] ;
}
DetJac = c11 * c22 - c12 * c21 ;
if (Element.Type == TRIANGLE)
Val_Surface += fabs(DetJac) * 0.5 ;
else if (Element.Type == QUADRANGLE)
Val_Surface += fabs(DetJac) * 4. ;
}
else {
Msg(GERROR, "Function 'SurfaceArea' only valid for Triangle or Quandrangle Elements");
}
}
}
Fct->Active->Case.SurfaceArea.Value = Val_Surface ;
}
V->Type = SCALAR ;
V->Val[0] = Fct->Active->Case.SurfaceArea.Value ;
V->Val[MAX_DIM] = 0;
for (k = 2 ; k < Current.NbrHar ; k += 2) {
V->Val[MAX_DIM* k] = V->Val[0] ;
V->Val[MAX_DIM* (k+1)] = 0 ;
}
GetDP_End ;
}
/* ------------------------------------------------------------------------ */
/* Transformation of a stiffness matrix (6x6) with 3 given angles */
/* ------------------------------------------------------------------------ */
#if defined(HAVE_GSL)
/* All the following should really be rewritten and generalized... */
void F_TransformTensor (F_ARG) {
Msg(GERROR, "Tuan's routines are not available with the GSL");
}
void F_TransformPerm (F_ARG) {
Msg(GERROR, "Tuan's routines are not available with the GSL");
}
void F_TransformPiezo (F_ARG) {
Msg(GERROR, "Tuan's routines are not available with the GSL");
}
void F_TransformPiezoT (F_ARG) {
Msg(GERROR, "Tuan's routines are not available with the GSL");
}
#else
void ludcmp1(double **a, int n, int *indx, double *d);
void lubksb1(double **a, int n, int *indx, double b[]);
/* Calculate the transformation stress tensor taking into account only one rotation
alpha/beta/gamma around the axis z/y/x respectively */
void T_Stress(double **T, int n, double alpha, double beta, double gamma)
{
int i,j,k;
double **T_z, **T_y, **T_x, **T_zy, sum;
T_z=dmatrix(0,n-1,0,n-1);
T_y=dmatrix(0,n-1,0,n-1);
T_x=dmatrix(0,n-1,0,n-1);
T_zy=dmatrix(0,n-1,0,n-1);
for (i=0;i<n;i++) {
for (j=0;j<n;j++) {
T_z[i][j] = 0.;
T_y[i][j] = 0.;
T_x[i][j] = 0.;
}
}
T_z[0][0] = SQU(cos(alpha));
T_z[0][1] = SQU(sin(alpha));
T_z[0][3] = 2*cos(alpha)*sin(alpha);
T_z[1][0] = SQU(sin(alpha));
T_z[1][1] = SQU(cos(alpha));
T_z[1][3] =-2*cos(alpha)*sin(alpha);
T_z[2][2] = 1.0;
T_z[3][0] =-cos(alpha)*sin(alpha);
T_z[3][1] = cos(alpha)*sin(alpha);
T_z[3][3] = SQU(cos(alpha))-SQU(sin(alpha));
T_z[4][4] = cos(alpha);
T_z[4][5] =-sin(alpha);
T_z[5][4] = sin(alpha);
T_z[5][5] = cos(alpha);
T_y[0][0] = SQU(cos(beta));
T_y[0][2] = SQU(sin(beta));
T_y[0][5] =-2*cos(beta)*sin(beta);
T_y[1][1] = 1.0;
T_y[2][0] = SQU(sin(beta));
T_y[2][2] = SQU(cos(beta));
T_y[2][5] = 2*cos(beta)*sin(beta);
T_y[3][3] = cos(beta);
T_y[3][4] =-sin(beta);
T_y[4][3] = sin(beta);
T_y[4][4] = cos(beta);
T_y[5][0] = cos(beta)*sin(beta);
T_y[5][2] =-cos(beta)*sin(beta);
T_y[5][5] = SQU(cos(beta))-SQU(sin(beta));
T_x[0][0] = 1.0;
T_x[1][1] = SQU(cos(gamma));
T_x[1][2] = SQU(sin(gamma));
T_x[1][4] = 2*cos(gamma)*sin(gamma);
T_x[2][1] = SQU(sin(gamma));
T_x[2][2] = SQU(cos(gamma));
T_x[2][4] =-2*cos(gamma)*sin(gamma);
T_x[3][3] = cos(gamma);
T_x[3][5] = sin(gamma);
T_x[4][1] =-cos(gamma)*sin(gamma);
T_x[4][2] = cos(gamma)*sin(gamma);
T_x[4][4] = SQU(cos(gamma))-SQU(sin(gamma));
T_x[5][3] =-sin(gamma);
T_x[5][5] = cos(gamma);
for (i=0;i<n;i++){
for (j=0;j<n;j++){
sum=0.0;
for(k=0;k<n;k++)
sum=sum+T_z[i][k]*T_y[k][j];
T_zy[i][j]=sum;
}
}
for (i=0;i<n;i++){
for (j=0;j<n;j++){
sum=0.0;
for(k=0;k<n;k++)
sum=sum+T_zy[i][k]*T_x[k][j];
T[i][j]=sum;
}
}
free_dmatrix(T_z, 0,n-1,0,n-1);
free_dmatrix(T_y, 0,n-1,0,n-1);
free_dmatrix(T_x, 0,n-1,0,n-1);
free_dmatrix(T_zy,0,n-1,0,n-1);
}
/* Calculate the transformation strain tensor taking into account only one rotation
alpha/beta/gamma around the axis z/y/x respectively */
void T_Strain(double **T, int n, double alpha, double beta, double gamma)
{
int i,j,k;
double **T_z, **T_y, **T_x, **T_zy, sum;
T_z=dmatrix(0,n-1,0,n-1);
T_y=dmatrix(0,n-1,0,n-1);
T_x=dmatrix(0,n-1,0,n-1);
T_zy=dmatrix(0,n-1,0,n-1);
for (i=0;i<n;i++) {
for (j=0;j<n;j++) {
T_z[i][j] = 0;
T_y[i][j] = 0;
T_x[i][j] = 0;
}
}
T_z[0][0] = SQU(cos(alpha));
T_z[0][1] = SQU(sin(alpha));
T_z[0][3] = cos(alpha)*sin(alpha);
T_z[1][0] = SQU(sin(alpha));
T_z[1][1] = SQU(cos(alpha));
T_z[1][3] =-cos(alpha)*sin(alpha);
T_z[2][2] = 1.0;
T_z[3][0] =-2*cos(alpha)*sin(alpha);
T_z[3][1] = 2*cos(alpha)*sin(alpha);
T_z[3][3] = SQU(cos(alpha))-SQU(sin(alpha));
T_z[4][4] = cos(alpha);
T_z[4][5] =-sin(alpha);
T_z[5][4] = sin(alpha);
T_z[5][5] = cos(alpha);
T_y[0][0] = SQU(cos(beta));
T_y[0][2] = SQU(sin(beta));
T_y[0][5] =-cos(beta)*sin(beta);
T_y[1][1] = 1.0;
T_y[2][0] = SQU(sin(beta));
T_y[2][2] = SQU(cos(beta));
T_y[2][5] = cos(beta)*sin(beta);
T_y[3][3] = cos(beta);
T_y[3][4] =-sin(beta);
T_y[4][3] = sin(beta);
T_y[4][4] = cos(beta);
T_y[5][0] = 2*cos(beta)*sin(beta);
T_y[5][2] =-2*cos(beta)*sin(beta);
T_y[5][5] = SQU(cos(beta))-SQU(sin(beta));
T_x[0][0] = 1.0;
T_x[1][1] = SQU(cos(gamma));
T_x[1][2] = SQU(sin(gamma));
T_x[1][4] = cos(gamma)*sin(gamma);
T_x[2][1] = SQU(sin(gamma));
T_x[2][2] = SQU(cos(gamma));
T_x[2][4] =-cos(gamma)*sin(gamma);
T_x[3][3] = cos(gamma);
T_x[3][5] = sin(gamma);
T_x[4][1] =-2*cos(gamma)*sin(gamma);
T_x[4][2] = 2*cos(gamma)*sin(gamma);
T_x[4][4] = SQU(cos(gamma))-SQU(sin(gamma));
T_x[5][3] =-sin(gamma);
T_x[5][5] = cos(gamma);
for (i=0;i<n;i++){
for (j=0;j<n;j++){
sum=0.0;
for(k=0;k<n;k++)
sum=sum+T_z[i][k]*T_y[k][j];
T_zy[i][j]=sum;
}
}
for (i=0;i<n;i++){
for (j=0;j<n;j++){
sum=0.0;
for(k=0;k<n;k++)
sum=sum+T_zy[i][k]*T_x[k][j];
T[i][j]=sum;
}
}
free_dmatrix(T_z, 0,n-1,0,n-1);
free_dmatrix(T_y, 0,n-1,0,n-1);
free_dmatrix(T_x, 0,n-1,0,n-1);
free_dmatrix(T_zy,0,n-1,0,n-1);
}
/* Tensor Transformation */
void F_TransformTensor (F_ARG) {
int ident, i, j, k, ii = 0, jj = 0, *indx, N = 6 ;
double **T_T, **T_S, **C, **C_trformed, **a, **y, d, *col, alpha, beta, Gamma;
GetDP_Begin("F_TransformTensor");
C = dmatrix(0,N-1,0,N-1);
C_trformed = dmatrix(0,N-1,0,N-1);
ident = (int)Fct->Para[0];
alpha = Fct->Para[1];
beta = Fct->Para[2];
Gamma = Fct->Para[3];
if( ( (A+0)->Type != TENSOR && (A+0)->Type != TENSOR_SYM && (A+0)->Type != TENSOR_DIAG ) ||
( (A+1)->Type != TENSOR && (A+1)->Type != TENSOR_SYM && (A+1)->Type != TENSOR_DIAG ) ||
( (A+2)->Type != TENSOR && (A+2)->Type != TENSOR_SYM && (A+2)->Type != TENSOR_DIAG ) ||
( (A+3)->Type != TENSOR && (A+3)->Type != TENSOR_SYM && (A+3)->Type != TENSOR_DIAG ) )
Msg(GERROR, "Function 'TransformTensor' requires 4 Tensors on input (NOT %s %s %s %s)",
Get_StringForDefine(Field_Type,(A+0)->Type),
Get_StringForDefine(Field_Type,(A+1)->Type),
Get_StringForDefine(Field_Type,(A+2)->Type),
Get_StringForDefine(Field_Type,(A+3)->Type) );
for(i=0;i<N;i++)
for(j=0;j<N;j++)
C[i][j] = 0.0; /* Reset C[i][j] to zeros */
for(i=0;i<4;i++) {
switch(i){
case 0 : ii=0; jj=0; break;
case 1 : ii=0; jj=3; break;
case 2 : ii=3; jj=0; break;
case 3 : ii=3; jj=3; break;
}
switch((A+i)->Type){
case TENSOR :
C[ii+0][jj+0] = (A+i)->Val[0]; C[ii+0][jj+1] = (A+i)->Val[1]; C[ii+0][jj+2] = (A+i)->Val[2];
C[ii+1][jj+0] = (A+i)->Val[3]; C[ii+1][jj+1] = (A+i)->Val[4]; C[ii+1][jj+2] = (A+i)->Val[5];
C[ii+2][jj+0] = (A+i)->Val[6]; C[ii+2][jj+1] = (A+i)->Val[7]; C[ii+2][jj+2] = (A+i)->Val[8];
break;
case TENSOR_DIAG :
C[ii+0][jj+0]=(A+i)->Val[0]; C[ii+1][jj+1]=(A+i)->Val[1]; C[ii+2][jj+2]=(A+i)->Val[2];
break;
case TENSOR_SYM :
C[ii+0][jj+0]=(A+i)->Val[0]; C[ii+0][jj+1]=(A+i)->Val[1]; C[ii+0][jj+2]=(A+i)->Val[2];
C[ii+1][jj+0]=C[ii+0][jj+1]; C[ii+1][jj+1]=(A+i)->Val[3]; C[ii+1][jj+2]=(A+i)->Val[4];
C[ii+2][jj+0]=C[ii+0][jj+2]; C[ii+2][jj+1]=C[ii+1][jj+2]; C[ii+2][jj+2]=(A+i)->Val[5];
break;
}
}
/* begin of transformation ! */
T_T=dmatrix(0,N-1,0,N-1);
T_S=dmatrix(0,N-1,0,N-1);
a=dmatrix(0,N-1,0,N-1);
y=dmatrix(0,N-1,0,N-1);
indx=ivector(0,N-1);
col=dvector(0,N-1);
T_Stress(T_T, N, alpha, beta, Gamma); /* Create the T_sigma */
T_Strain(T_S, N, alpha, beta, Gamma); /* Create the T_epsilon */
ludcmp1(T_T,N,indx,&d); /* LU decomposition of T_sigma */
for(j=0;j<N;j++) {
for(i=0;i<N;i++) col[i]=0.0;
col[j]=1.0;
lubksb1(T_T,N,indx,col); /* Forward and back substitution */
for(i=0;i<N;i++) y[i][j]= col[i]; /* inverse( T_sigma ) = y */
}
for (i=0;i<N;i++){
for (j=0;j<N;j++){
a[i][j]=0.0;
for(k=0;k<N;k++) a[i][j]=a[i][j]+y[i][k]*C[k][j]; /* inverse(T_sigma)*C[i][j] */
}
}
for (i=0;i<N;i++){
for (j=0;j<N;j++){
C_trformed[i][j]=0.0;
for(k=0;k<N;k++) C_trformed[i][j]=C_trformed[i][j]+a[i][k]*T_S[k][j];
/* C[i][j]' = inverse(T_sigma)*C[i][j]*T_epsilon */
}
}
switch (ident) {
case 11 :
(V+0)->Val[0] = C_trformed[0][0]; (V+0)->Val[1] = C_trformed[0][1]; (V+0)->Val[2] = C_trformed[0][2];
(V+0)->Val[3] = C_trformed[1][0]; (V+0)->Val[4] = C_trformed[1][1]; (V+0)->Val[5] = C_trformed[1][2];
(V+0)->Val[6] = C_trformed[2][0]; (V+0)->Val[7] = C_trformed[2][1]; (V+0)->Val[8] = C_trformed[2][2];
/*
fprintf(stderr,"(C11 0) %.16g %.16g %.16g\n", C_trformed[0][0], C_trformed[0][1], C_trformed[0][2]);
fprintf(stderr,"(C11 1) %.16g %.16g %.16g\n", C_trformed[1][0], C_trformed[1][1], C_trformed[1][2]);
fprintf(stderr,"(C11 2) %.16g %.16g %.16g\n", C_trformed[2][0], C_trformed[2][1], C_trformed[2][2]);
*/
break ;
case 12 :
(V+0)->Val[0] = C_trformed[0][3]; (V+0)->Val[1] = C_trformed[0][4]; (V+0)->Val[2] = C_trformed[0][5];
(V+0)->Val[3] = C_trformed[1][3]; (V+0)->Val[4] = C_trformed[1][4]; (V+0)->Val[5] = C_trformed[1][5];
(V+0)->Val[6] = C_trformed[2][3]; (V+0)->Val[7] = C_trformed[2][4]; (V+0)->Val[8] = C_trformed[2][5];
/*
fprintf(stderr,"(C12 0) %.16g %.16g %.16g\n", C_trformed[0][3], C_trformed[0][4], C_trformed[0][5]);
fprintf(stderr,"(C12 1) %.16g %.16g %.16g\n", C_trformed[1][3], C_trformed[1][4], C_trformed[1][5]);
fprintf(stderr,"(C12 2) %.16g %.16g %.16g\n", C_trformed[2][3], C_trformed[2][4], C_trformed[2][5]);
*/
break ;
case 21 :
(V+0)->Val[0] = C_trformed[3][0]; (V+0)->Val[1] = C_trformed[3][1]; (V+0)->Val[2] = C_trformed[3][2];
(V+0)->Val[3] = C_trformed[4][0]; (V+0)->Val[4] = C_trformed[4][1]; (V+0)->Val[5] = C_trformed[4][2];
(V+0)->Val[6] = C_trformed[5][0]; (V+0)->Val[7] = C_trformed[5][1]; (V+0)->Val[8] = C_trformed[5][2];
/*
fprintf(stderr,"(C21 0) %.16g %.16g %.16g\n", C_trformed[3][0], C_trformed[3][1], C_trformed[3][2]);
fprintf(stderr,"(C21 1) %.16g %.16g %.16g\n", C_trformed[4][0], C_trformed[4][1], C_trformed[4][2]);
fprintf(stderr,"(C21 2) %.16g %.16g %.16g\n", C_trformed[5][0], C_trformed[5][1], C_trformed[5][2]);
*/
break ;
case 22 :
(V+0)->Val[0] = C_trformed[3][3]; (V+0)->Val[1] = C_trformed[3][4]; (V+0)->Val[2] = C_trformed[3][5];
(V+0)->Val[3] = C_trformed[4][3]; (V+0)->Val[4] = C_trformed[4][4]; (V+0)->Val[5] = C_trformed[4][5];
(V+0)->Val[6] = C_trformed[5][3]; (V+0)->Val[7] = C_trformed[5][4]; (V+0)->Val[8] = C_trformed[5][5];
/*
fprintf(stderr,"(C22 0) %.16g %.16g %.16g\n", C_trformed[3][3], C_trformed[3][4], C_trformed[3][5]);
fprintf(stderr,"(C22 1) %.16g %.16g %.16g\n", C_trformed[4][3], C_trformed[4][4], C_trformed[4][5]);
fprintf(stderr,"(C22 2) %.16g %.16g %.16g\n", C_trformed[5][3], C_trformed[5][4], C_trformed[5][5]);
*/
break ;
}
free_dmatrix(C, 0,N-1,0,N-1);
free_dmatrix(C_trformed, 0,N-1,0,N-1);
free_dmatrix(T_T, 0,N-1,0,N-1);
free_dmatrix(T_S, 0,N-1,0,N-1);
free_dmatrix(a, 0,N-1,0,N-1);
free_dmatrix(y, 0,N-1,0,N-1);
free_ivector(indx, 0,N-1);
free_dvector(col, 0,N-1);
V->Type = TENSOR; /* We do not know exactly the type of tensor,
so TENSOR is a most general choice */
/* end of transformation */
GetDP_End ;
}
/* Calculate the transformation permitivity tensor taking into account only one rotation
alpha/beta/gamma around the axis z/y/x respectively */
void T_Epsilon(double **T, double alpha, double beta, double gamma)
{
T[0][0] = cos(alpha)* cos(beta);
T[0][1] = sin(alpha)*cos(gamma)+cos(alpha)*sin(beta)*sin(gamma);
T[0][2] = sin(alpha)*sin(gamma)-cos(alpha)*sin(beta)*cos(gamma);
T[1][0] =-sin(alpha)*cos(beta);
T[1][1] = cos(alpha)*cos(gamma)-sin(alpha)*sin(beta)*sin(gamma);
T[1][2] = cos(alpha)*sin(gamma)+sin(alpha)*sin(beta)*cos(gamma);
T[2][0] = sin(beta);
T[2][1] =-cos(beta)*sin(gamma);
T[2][2] = cos(beta)*cos(gamma);
}
void F_TransformPerm (F_ARG) {
int i, j, k, N=3 ;
double **T_eps,**EPSr,**EPSr_trformed,**a,**y,alpha,beta,Gamma;
GetDP_Begin("F_TransformPerm");
EPSr = dmatrix(0,N-1,0,N-1);
EPSr_trformed = dmatrix(0,N-1,0,N-1);
alpha = Fct->Para[0];
beta = Fct->Para[1];
Gamma = Fct->Para[2];
if ( A->Type != TENSOR && A->Type != TENSOR_SYM && A->Type != TENSOR_DIAG )
Msg(GERROR, "Wrong type of argument for function 'TransformTensor2' (NOT %s) ",
Get_StringForDefine(Field_Type,A->Type));
for(i=0;i<N;i++)
for(j=0;j<N;j++)
EPSr[i][j] = 0.0; /* Reset EPSr[i][j] to zeros */
switch(A->Type){
case TENSOR :
EPSr[0][0] = A->Val[0]; EPSr[0][1] = A->Val[1]; EPSr[0][2] = A->Val[2];
EPSr[1][0] = A->Val[3]; EPSr[1][1] = A->Val[4]; EPSr[1][2] = A->Val[5];
EPSr[2][0] = A->Val[6]; EPSr[2][1] = A->Val[7]; EPSr[2][2] = A->Val[8];
break;
case TENSOR_DIAG :
EPSr[0][0] = A->Val[0]; EPSr[1][1] = A->Val[1]; EPSr[2][2]=A->Val[2];
break;
case TENSOR_SYM :
EPSr[0][0] = A->Val[0]; EPSr[0][1] = A->Val[1]; EPSr[0][2]=A->Val[2];
EPSr[1][0] = EPSr[0][1]; EPSr[1][1] = A->Val[3]; EPSr[1][2]=A->Val[4];
EPSr[2][0] = EPSr[0][2]; EPSr[2][1] = EPSr[1][2]; EPSr[2][2]=A->Val[5];
break;
}
/* begin of transformation ! */
T_eps=dmatrix(0,N-1,0,N-1);
a=dmatrix(0,N-1,0,N-1);
y=dmatrix(0,N-1,0,N-1);
T_Epsilon(T_eps, alpha, beta, Gamma);
for(i=0;i<N;i++){
for(j=0;j<N;j++){
if (i!=j) y[i][j]=T_eps[j][i];
else y[i][j]=T_eps[i][j];
}
}
for(i=0;i<N;i++){
for (j=0;j<N;j++){
a[i][j]=0.0;
for(k=0;k<N;k++)
a[i][j]=a[i][j]+y[i][k]*EPSr[k][j];
}
}
for(i=0;i<N;i++){
for(j=0;j<N;j++){
EPSr_trformed[i][j]=0.0;
for(k=0;k<N;k++)
EPSr_trformed[i][j]=EPSr_trformed[i][j]+a[i][k]*T_eps[k][j];
}
}
V->Val[0] = EPSr_trformed[0][0]; V->Val[1] = EPSr_trformed[0][1]; V->Val[2] = EPSr_trformed[0][2];
V->Val[3] = EPSr_trformed[1][0]; V->Val[4] = EPSr_trformed[1][1]; V->Val[5] = EPSr_trformed[1][2];
V->Val[6] = EPSr_trformed[2][0]; V->Val[7] = EPSr_trformed[2][1]; V->Val[8] = EPSr_trformed[2][2];
/*
fprintf(stderr,"(epsr 0) %.16g %.16g %.16g\n", EPSr_trformed[0][0], EPSr_trformed[0][1], EPSr_trformed[0][2]);
fprintf(stderr,"(epsr 1) %.16g %.16g %.16g\n", EPSr_trformed[1][0], EPSr_trformed[1][1], EPSr_trformed[1][2]);
fprintf(stderr,"(epsr 2) %.16g %.16g %.16g\n", EPSr_trformed[2][0], EPSr_trformed[2][1], EPSr_trformed[2][2]);
*/
free_dmatrix(EPSr, 0,N-1,0,N-1);
free_dmatrix(EPSr_trformed, 0,N-1,0,N-1);
free_dmatrix(T_eps, 0,N-1,0,N-1);
free_dmatrix(a, 0,N-1,0,N-1);
free_dmatrix(y, 0,N-1,0,N-1);
V->Type = TENSOR;
/* end of transformation */
GetDP_End ;
}
void F_TransformPiezo (F_ARG) {
int ident, i, j, k, ii = 0, jj = 0, N = 6 ;
double **T_eps, **T_S, **e, **e_trformed, **a, **y, alpha, beta, Gamma;
GetDP_Begin("F_TransformPiezo");
e = dmatrix(0,N/2-1,0,N-1);
e_trformed = dmatrix(0,N/2-1,0,N-1);
ident = (int)Fct->Para[0];
alpha = Fct->Para[1];
beta = Fct->Para[2];
Gamma = Fct->Para[3];
if( ( (A+0)->Type != TENSOR && (A+0)->Type != TENSOR_SYM && (A+0)->Type != TENSOR_DIAG ) ||
( (A+1)->Type != TENSOR && (A+1)->Type != TENSOR_SYM && (A+1)->Type != TENSOR_DIAG ) )
Msg(GERROR, "Function 'TransformTensor' requires 2 Tensors on input (NOT %s %s )",
Get_StringForDefine(Field_Type,(A+0)->Type),
Get_StringForDefine(Field_Type,(A+1)->Type) );
for(i=0;i<N/2;i++)
for(j=0;j<N;j++)
e[i][j] = 0.0; /* Reset e[i][j] to zeros */
for(i=0;i<2;i++) {
switch(i){
case 0 : ii=0; jj=0; break;
case 1 : ii=0; jj=3; break;
}
switch((A+i)->Type){
case TENSOR :
e[ii+0][jj+0] = (A+i)->Val[0]; e[ii+0][jj+1] = (A+i)->Val[1]; e[ii+0][jj+2] = (A+i)->Val[2];
e[ii+1][jj+0] = (A+i)->Val[3]; e[ii+1][jj+1] = (A+i)->Val[4]; e[ii+1][jj+2] = (A+i)->Val[5];
e[ii+2][jj+0] = (A+i)->Val[6]; e[ii+2][jj+1] = (A+i)->Val[7]; e[ii+2][jj+2] = (A+i)->Val[8];
break;
case TENSOR_DIAG :
e[ii+0][jj+0]=(A+i)->Val[0]; e[ii+1][jj+1]=(A+i)->Val[1]; e[ii+2][jj+2]=(A+i)->Val[2];
break;
case TENSOR_SYM :
e[ii+0][jj+0]=(A+i)->Val[0]; e[ii+0][jj+1]=(A+i)->Val[1]; e[ii+0][jj+2]=(A+i)->Val[2];
e[ii+1][jj+0]=e[ii+0][jj+1]; e[ii+1][jj+1]=(A+i)->Val[3]; e[ii+1][jj+2]=(A+i)->Val[4];
e[ii+2][jj+0]=e[ii+0][jj+2]; e[ii+2][jj+1]=e[ii+1][jj+2]; e[ii+2][jj+2]=(A+i)->Val[5];
break;
}
}
/* begin of transformation ! */
T_S=dmatrix(0,N-1,0,N-1);
T_eps=dmatrix(0,N/2-1,0,N/2-1);
a=dmatrix(0,N/2-1,0,N-1);
y=dmatrix(0,N/2-1,0,N/2-1);
T_Strain(T_S, N, alpha, beta, Gamma);
T_Epsilon(T_eps, alpha, beta, Gamma);
for(i=0;i<N/2;i++){
for(j=0;j<N/2;j++){
if (i!=j) y[i][j]=T_eps[j][i];
else y[i][j]=T_eps[i][j];
}
}
for(i=0;i<N/2;i++){
for (j=0;j<N;j++){
a[i][j]=0.0;
for(k=0;k<N/2;k++)
a[i][j]=a[i][j]+y[i][k]*e[k][j];
}
}
for(i=0;i<N/2;i++){
for(j=0;j<N;j++){
e_trformed[i][j]=0.0;
for(k=0;k<N;k++)
e_trformed[i][j]=e_trformed[i][j]+a[i][k]*T_S[k][j];
}
}
switch (ident) {
case 1 :
V->Val[0] = e_trformed[0][0]; V->Val[1] = e_trformed[0][1]; V->Val[2] = e_trformed[0][2];
V->Val[3] = e_trformed[1][0]; V->Val[4] = e_trformed[1][1]; V->Val[5] = e_trformed[1][2];
V->Val[6] = e_trformed[2][0]; V->Val[7] = e_trformed[2][1]; V->Val[8] = e_trformed[2][2];
/*
fprintf(stderr,"(e 11 0) %.16g %.16g %.16g\n", e_trformed[0][0], e_trformed[0][1], e_trformed[0][2]);
fprintf(stderr,"(e 11 1) %.16g %.16g %.16g\n", e_trformed[1][0], e_trformed[1][1], e_trformed[1][2]);
fprintf(stderr,"(e 11 2) %.16g %.16g %.16g\n", e_trformed[2][0], e_trformed[2][1], e_trformed[2][2]);
*/
break ;
case 2 :
V->Val[0] = e_trformed[0][3]; V->Val[1] = e_trformed[0][4]; V->Val[2] = e_trformed[0][5];
V->Val[3] = e_trformed[1][3]; V->Val[4] = e_trformed[1][4]; V->Val[5] = e_trformed[1][5];
V->Val[6] = e_trformed[2][3]; V->Val[7] = e_trformed[2][4]; V->Val[8] = e_trformed[2][5];
/*
fprintf(stderr,"(e 12 0) %.16g %.16g %.16g\n", e_trformed[0][3], e_trformed[0][4], e_trformed[0][5]);
fprintf(stderr,"(e 12 1) %.16g %.16g %.16g\n", e_trformed[1][3], e_trformed[1][4], e_trformed[1][5]);
fprintf(stderr,"(e 12 2) %.16g %.16g %.16g\n", e_trformed[2][3], e_trformed[2][4], e_trformed[2][5]);
*/
break ;
}
free_dmatrix(e, 0,N/2-1,0,N-1);
free_dmatrix(e_trformed, 0,N/2-1,0,N-1);
free_dmatrix(T_S, 0,N-1,0,N-1);
free_dmatrix(T_eps, 0,N/2-1,0,N/2-1);
free_dmatrix(a, 0,N/2-1,0,N-1);
free_dmatrix(y, 0,N/2-1,0,N/2-1);
V->Type = TENSOR;
/* end of transformation */
GetDP_End ;
}
void F_TransformPiezoT (F_ARG) {
int ident, i, j, k,*indx, ii = 0, jj = 0, N=6 ;
double **T_eps, **T_T, **e, **eT, **eT_trformed, d, *col, **a, **y, alpha, beta, Gamma;
GetDP_Begin("F_TransformCij");
e = dmatrix(0,N/2-1,0,N-1);
eT = dmatrix(0,N-1,0,N/2-1);
eT_trformed = dmatrix(0,N-1,0,N/2-1);
ident = (int)Fct->Para[0];
alpha = Fct->Para[1];
beta = Fct->Para[2];
Gamma = Fct->Para[3];
if( ( (A+0)->Type != TENSOR && (A+0)->Type != TENSOR_SYM && (A+0)->Type != TENSOR_DIAG ) ||
( (A+1)->Type != TENSOR && (A+3)->Type != TENSOR_SYM && (A+3)->Type != TENSOR_DIAG ) )
Msg(GERROR, "Function 'TransformTensor' requires 2 Tensors on input (NOT %s %s )",
Get_StringForDefine(Field_Type,(A+0)->Type),
Get_StringForDefine(Field_Type,(A+1)->Type) );
for(i=0;i<N/2;i++)
for(j=0;j<N;j++)
e[i][j] = 0.0; /* Reset e[i][j] to zeros */
for(i=0;i<2;i++) {
switch(i){
case 0 : ii=0; jj=0; break;
case 1 : ii=0; jj=3; break;
}
switch((A+i)->Type){
case TENSOR :
e[ii+0][jj+0] = (A+i)->Val[0]; e[ii+0][jj+1] = (A+i)->Val[1]; e[ii+0][jj+2] = (A+i)->Val[2];
e[ii+1][jj+0] = (A+i)->Val[3]; e[ii+1][jj+1] = (A+i)->Val[4]; e[ii+1][jj+2] = (A+i)->Val[5];
e[ii+2][jj+0] = (A+i)->Val[6]; e[ii+2][jj+1] = (A+i)->Val[7]; e[ii+2][jj+2] = (A+i)->Val[8];
break;
case TENSOR_DIAG :
e[ii+0][jj+0]=(A+i)->Val[0]; e[ii+1][jj+1]=(A+i)->Val[1]; e[ii+2][jj+2]=(A+i)->Val[2];
break;
case TENSOR_SYM :
e[ii+0][jj+0]=(A+i)->Val[0]; e[ii+0][jj+1]=(A+i)->Val[1]; e[ii+0][jj+2]=(A+i)->Val[2];
e[ii+1][jj+0]=e[ii+0][jj+1]; e[ii+1][jj+1]=(A+i)->Val[3]; e[ii+1][jj+2]=(A+i)->Val[4];
e[ii+2][jj+0]=e[ii+0][jj+2]; e[ii+2][jj+1]=e[ii+1][jj+2]; e[ii+2][jj+2]=(A+i)->Val[5];
break;
}
}
/* begin of transformation ! */
T_T=dmatrix(0,N-1,0,N-1);
T_eps=dmatrix(0,N/2-1,0,N/2-1);
a=dmatrix(0,N-1,0,N/2-1);
y=dmatrix(0,N-1,0,N-1);
indx=ivector(0,N-1);
col=dvector(0,N-1);
T_Stress(T_T, N, alpha, beta, Gamma);
T_Epsilon(T_eps, alpha, beta, Gamma);
ludcmp1(T_T,N,indx,&d); /* LU decomposition of T_sigma */
for(j=0;j<N;j++) {
for(i=0;i<N;i++) col[i]=0.0;
col[j]=1.0;
lubksb1(T_T,N,indx,col); /* Forward and back substitution */
for(i=0;i<N;i++) y[i][j]= col[i]; /* inverse( T_sigma ) = y */
}
for(i=0;i<N;i++){
for(j=0;j<N/2;j++){
if (i!=j) eT[i][j]=e[j][i];
else eT[i][j]=e[i][j];
}
}
for (i=0;i<N;i++){
for (j=0;j<N/2;j++){
a[i][j]=0.0;
for(k=0;k<N;k++) a[i][j]=a[i][j]+y[i][k]*eT[k][j];
}
}
for (i=0;i<N;i++){
for (j=0;j<N/2;j++){
eT_trformed[i][j]=0.0;
for(k=0;k<N/2;k++) eT_trformed[i][j]=eT_trformed[i][j]+a[i][k]*T_eps[k][j];
}
}
switch (ident) {
case 1:
V->Val[0] = eT_trformed[0][0]; V->Val[1] = eT_trformed[0][1]; V->Val[2] = eT_trformed[0][2];
V->Val[3] = eT_trformed[1][0]; V->Val[4] = eT_trformed[1][1]; V->Val[5] = eT_trformed[1][2];
V->Val[6] = eT_trformed[2][0]; V->Val[7] = eT_trformed[2][1]; V->Val[8] = eT_trformed[2][2];
/*
fprintf(stderr,"(eT 11 0) %.16g %.16g %.16g\n", eT_trformed[0][0], eT_trformed[0][1], eT_trformed[0][2]);
fprintf(stderr,"(eT 11 1) %.16g %.16g %.16g\n", eT_trformed[1][0], eT_trformed[1][1], eT_trformed[1][2]);
fprintf(stderr,"(eT 11 2) %.16g %.16g %.16g\n", eT_trformed[2][0], eT_trformed[2][1], eT_trformed[2][2]);
*/
break ;
case 2 :
V->Val[0] = eT_trformed[3][0]; V->Val[1] = eT_trformed[3][1]; V->Val[2] = eT_trformed[3][2];
V->Val[3] = eT_trformed[4][0]; V->Val[4] = eT_trformed[4][1]; V->Val[5] = eT_trformed[4][2];
V->Val[6] = eT_trformed[5][0]; V->Val[7] = eT_trformed[5][1]; V->Val[8] = eT_trformed[5][2];
/*
fprintf(stderr,"(eT 21 0) %.16g %.16g %.16g\n", eT_trformed[3][0], eT_trformed[3][1], eT_trformed[3][2]);
fprintf(stderr,"(eT 21 1) %.16g %.16g %.16g\n", eT_trformed[4][0], eT_trformed[4][1], eT_trformed[4][2]);
fprintf(stderr,"(eT 21 2) %.16g %.16g %.16g\n", eT_trformed[5][0], eT_trformed[5][1], eT_trformed[5][2]);
*/
break ;
}
free_dmatrix(e, 0,N/2-1,0,N-1);
free_dmatrix(eT, 0,N/2,0,N/2-1);
free_dmatrix(eT_trformed, 0,N-1,0,N/2-1);
free_dmatrix(T_T, 0,N-1,0,N-1);
free_dmatrix(T_eps, 0,N/2-1,0,N/2-1);
free_dmatrix(a, 0,N-1,0,N/2-1);
free_dmatrix(y, 0,N-1,0,N-1);
free_ivector(indx, 0,N-1);
free_dvector(col, 0,N-1);
V->Type = TENSOR;
/* end of transformation */
GetDP_End ;
}
#endif /* #if !defined(HAVE_GSL) */
/* ------------------------------------------------------------------------ */
/* Start: V I R T U A L W O R K */
/* ------------------------------------------------------------------------ */
void JacobianVol_dx (struct Element * Element,
MATRIX3x3 * Jac, double DetJac,
MATRIX3x3 * Jac_dx, double * DetJac_dx, int i) {
GetDP_Begin("JacobianVol2D_dx");
Jac_dx->c11 = Element->dndu[i][0] ; Jac_dx->c12 = Element->dndu[i][0] ; Jac_dx->c13 = Element->dndu[i][0] ;
Jac_dx->c21 = Element->dndu[i][1] ; Jac_dx->c22 = Element->dndu[i][1] ; Jac_dx->c23 = Element->dndu[i][1] ;
Jac_dx->c31 = Element->dndu[i][2] ; Jac_dx->c32 = Element->dndu[i][2] ; Jac_dx->c33 = Element->dndu[i][2] ;
DetJac_dx[0] = Jac_dx->c11 * ( Jac->c22 * Jac->c33 - Jac->c23 * Jac->c32 )
- Jac_dx->c21 * ( Jac->c12 * Jac->c33 - Jac->c13 * Jac->c32 )
+ Jac_dx->c31 * ( Jac->c12 * Jac->c23 - Jac->c22 * Jac->c13 );
DetJac_dx[1] = - Jac_dx->c12 * ( Jac->c21 * Jac->c33 - Jac->c23 * Jac->c31 )
+ Jac_dx->c22 * ( Jac->c11 * Jac->c33 - Jac->c13 * Jac->c31 )
- Jac_dx->c32 * ( Jac->c11 * Jac->c23 - Jac->c13 * Jac->c21 );
DetJac_dx[2] = Jac_dx->c11 * ( Jac->c21 * Jac->c32 - Jac->c22 * Jac->c31 )
- Jac_dx->c21 * ( Jac->c11 * Jac->c32 - Jac->c12 * Jac->c31 )
+ Jac_dx->c31 * ( Jac->c11 * Jac->c22 - Jac->c12 * Jac->c21 );
/* DetJac_dx[0] = Jac_dx->c11 * Jac->c22 - Jac_dx->c21 * Jac->c12 ; */
/* DetJac_dx[1] = Jac_dx->c22 * Jac->c11 - Jac_dx->c12 * Jac->c21 ; */
/* DetJac_dx[2] = 0. ; */
if (DetJac < 0) {
DetJac_dx[0] *= -1.;
DetJac_dx[1] *= -1.;
DetJac_dx[2] *= -1.;
}
GetDP_End ;
}
/* F_VirtualWork */
void F_VirtualWork (F_ARG) {
int i, numNode, Type_Dimension;
double u, v, w;
MATRIX3x3 Jac;
MATRIX3x3 Jac_dx;
double DetJac;
double DetJac_dx [3];
double valField[3], s[3] ;
GetDP_Begin("F_VirtualWork");
/* numNode = (int)((A+1)->Val[0]); */
numNode = Current.NumEntity;
i = 0;
while (i < Current.Element->GeoElement->NbrNodes &&
Current.Element->GeoElement->NumNodes[i]!=numNode) i++;
if (i < Current.Element->GeoElement->NbrNodes ) {
/*
printf("element : %d %d\n", Current.Element->GeoElement->Num, i+1);
*/
valField[0] = A->Val[0];
valField[1] = A->Val[1];
valField[2] = A->Val[2];
u = Current.u; v = Current.v; w = Current.w;
Get_BFGeoElement(Current.Element, u,v,w);
Type_Dimension = Current.GeoData->Dimension;
switch (Type_Dimension) {
case _2D :
DetJac = JacobianVol2D (Current.Element, &Jac) ;
JacobianVol_dx (Current.Element, &Jac, DetJac, &Jac_dx, DetJac_dx, i) ;
s[0] = DetJac_dx[0] * ( - valField[0] * valField[0] +
valField[1] * valField[1] +
valField[2] * valField[2] )
- 2 * DetJac_dx[1] * valField[0] * valField[1]
- 2 * DetJac_dx[2] * valField[0] * valField[2];
/* DetJac_dx[0] * (valField[1] * valField[1] - valField[0] * valField[0]) */
/* - 2 * DetJac_dx[1] * valField[0] * valField[1] ; */
s[1] = DetJac_dx[1] * ( valField[0] * valField[0] -
valField[1] * valField[1] +
valField[2] * valField[2] )
- 2 * DetJac_dx[0] * valField[1] * valField[0]
- 2 * DetJac_dx[2] * valField[1] * valField[2];
/* DetJac_dx[1] * (valField[0] * valField[0] - valField[1] * valField[1]) */
/* - 2 * DetJac_dx[0] * valField[0] * valField[1] ; */
s[2] = DetJac_dx[2] * ( valField[0] * valField[0] +
valField[1] * valField[1] -
valField[2] * valField[2] )
- 2 * DetJac_dx[0] * valField[2] * valField[0]
- 2 * DetJac_dx[1] * valField[2] * valField[1];
s[0] /= fabs(DetJac);
s[1] /= fabs(DetJac);
s[2] /= fabs(DetJac);
break;
case _3D :
DetJac = JacobianVol3D (Current.Element, &Jac) ;
JacobianVol_dx (Current.Element, &Jac, DetJac, &Jac_dx, DetJac_dx, i) ;
s[0] = DetJac_dx[0] * ( - valField[0] * valField[0] +
valField[1] * valField[1] +
valField[2] * valField[2] )
- 2 * DetJac_dx[1] * valField[0] * valField[1]
- 2 * DetJac_dx[2] * valField[0] * valField[2];
/* DetJac_dx[0] * (valField[1] * valField[1] - valField[0] * valField[0]) */
/* - 2 * DetJac_dx[1] * valField[0] * valField[1] ; */
s[1] = DetJac_dx[1] * ( valField[0] * valField[0] -
valField[1] * valField[1] +
valField[2] * valField[2] )
- 2 * DetJac_dx[0] * valField[1] * valField[0]
- 2 * DetJac_dx[2] * valField[1] * valField[2];
/* DetJac_dx[1] * (valField[0] * valField[0] - valField[1] * valField[1]) */
/* - 2 * DetJac_dx[0] * valField[0] * valField[1] ; */
s[2] = DetJac_dx[2] * ( valField[0] * valField[0] +
valField[1] * valField[1] -
valField[2] * valField[2] )
- 2 * DetJac_dx[0] * valField[2] * valField[0]
- 2 * DetJac_dx[1] * valField[2] * valField[1];
s[0] /= fabs(DetJac);
s[1] /= fabs(DetJac);
s[2] /= fabs(DetJac);
break;
default :
s[0] = 0.; s[1] = 0.; s[2] = 0.;
break;
}
}
else {
s[0] = 0.; s[1] = 0.; s[2] = 0.;
}
V->Type = VECTOR ;
V->Val[0] = s[0] ;
V->Val[1] = s[1] ;
V->Val[2] = s[2] ;
GetDP_End ;
}
#undef F_ARG
syntax highlighted by Code2HTML, v. 0.9.1