v0.6.10
Macros | Functions | Variables
quality_tetrahedorn.c File Reference

Go to the source code of this file.

Macros

#define QUALITY_VOLUME_LENGTH
 

Functions

PetscErrorCode get_edges_from_elem_coords (double *coords, double *coords_edges)
 
PetscErrorCode calculate_lrms (double *dofs_egdes_X, double *lrms)
 
PetscErrorCode calculate_push_edge_relative (double *edge_coords, __CLPK_doublecomplex *xH, __CLPK_doublecomplex *dofs_egdes_X, __CLPK_doublecomplex *xlrms)
 
PetscErrorCode calculate_xQ (__CLPK_doublecomplex *xlrms, __CLPK_doublecomplex *inv_xH, __CLPK_doublecomplex *dofs_egdes_X_1, __CLPK_doublecomplex *dofs_egdes_X_2, __CLPK_doublecomplex *xQ)
 
void set_qual_ver (int ver)
 
int get_qual_ver ()
 
PetscErrorCode quality_volume_length_F (double V, double *alpha2, double gamma, double *diffN, double *coords_edges, double *dofs_X, double *dofs_x, double *dofs_iX, double *dofs_ix, double *quality0, double *quality, double *b, double *F, double *iF)
 
int quality_volume_length_K (double eps, double V, double *alpha2, double gamma, double *diffN, double *coords_edges, double *dofs_X, double *dofs_x, double *K, double *Koff)
 

Variables

static PetscErrorCode ierr
 
static int qual_ver = 1
 

Macro Definition Documentation

◆ QUALITY_VOLUME_LENGTH

#define QUALITY_VOLUME_LENGTH
Value:
__CLPK_doublecomplex dofs_egdes_X[3*6]; \
ierr = calculate_push_edge_relative(coords_edges,xH,dofs_egdes_X,&xlrms); CHKERRQ(ierr); \
/* calculeate xQ */ \
__CLPK_doublecomplex xQ[9]; \
ierr = calculate_xQ(&xlrms,inv_xH,dofs_egdes_Chi,dofs_egdes_X,xQ); CHKERRQ(ierr); \
/* some useful varibles */ \
__CLPK_doublecomplex xV; \
xV.r = det_xH.r*V; \
xV.i = det_xH.i*V; \
complex double complex_b = (det_xH.r+I*det_xH.i)/cpow((xlrms.r+I*xlrms.i)/xlrms0.r,3); \
__CLPK_doublecomplex xb = { creal(complex_b), cimag(complex_b) }; \
complex double complex_q = 6.*sqrt(2.)*(xV.r+I*xV.i)/cpow(xlrms.r+I*xlrms.i,3.); \
__CLPK_doublecomplex xq = { creal(complex_q), cimag(complex_q) }; \
complex double complex_grad; \
if( qual_ver == 0 ) { \
/* barrier and quality gradient change */ \
complex_grad = (xb.r+I*xb.i)/(1.-gamma)-1./((xb.r+I*xb.i)-gamma); \
} \
if( qual_ver == 1 ) { \
/* barrier and quality gradient */ \
complex_grad = (xq.r+I*xq.i)/(1-gamma)-1./((xq.r+I*xq.i)-gamma); \
} \
if( qual_ver == 2 ) { \
/* quality gradient */ \
complex_grad = xq.r+I*xq.i; \
} \
if( qual_ver == 3 ) { \
/* quality gradient scaled by volume */ \
complex_grad = (xV.r+I*xV.i)*((xb.r+I*xb.i)/(1.-gamma)-1./((xb.r+I*xb.i)-gamma)); \
} \
__CLPK_doublecomplex xgrad = { creal(complex_grad), cimag(complex_grad) }; \
cblas_zscal(9,&xgrad,xQ,1);
static PetscErrorCode ierr
static int qual_ver
PetscErrorCode calculate_xQ(__CLPK_doublecomplex *xlrms, __CLPK_doublecomplex *inv_xH, __CLPK_doublecomplex *dofs_egdes_X_1, __CLPK_doublecomplex *dofs_egdes_X_2, __CLPK_doublecomplex *xQ)
__CLPK_doublereal i
Definition: lapack_wrap.h:45
CHKERRQ(ierr)
__CLPK_doublereal r
Definition: lapack_wrap.h:45
PetscErrorCode calculate_push_edge_relative(double *edge_coords, __CLPK_doublecomplex *xH, __CLPK_doublecomplex *dofs_egdes_X, __CLPK_doublecomplex *xlrms)

Definition at line 119 of file quality_tetrahedorn.c.

Function Documentation

◆ calculate_lrms()

PetscErrorCode calculate_lrms ( double dofs_egdes_X,
double lrms 
)

Definition at line 33 of file quality_tetrahedorn.c.

33  {
34  PetscFunctionBegin;
35  *lrms = 0;
36  int ee = 0;
37  //loop over edges
38  for(;ee<6;ee++) {
39  double edge_relative[3];
40  bzero(edge_relative,3*sizeof(double));
41  int jj = 0;
42  for(;jj<3;jj++) {
43  edge_relative[jj] = dofs_egdes_X[6*ee+jj];
44  edge_relative[jj] -= dofs_egdes_X[6*ee+3+jj];
45  *lrms += pow(edge_relative[jj],2.);
46  }
47  }
48  *lrms = sqrt( (1./6.)*(*lrms) );
49  PetscFunctionReturn(0);
50 }

◆ calculate_push_edge_relative()

PetscErrorCode calculate_push_edge_relative ( double edge_coords,
__CLPK_doublecomplex xH,
__CLPK_doublecomplex dofs_egdes_X,
__CLPK_doublecomplex xlrms 
)

Definition at line 51 of file quality_tetrahedorn.c.

51  {
52  PetscFunctionBegin;
53  bzero(dofs_egdes_X,3*6*sizeof(__CLPK_doublecomplex));
54  double complex lrms = 0;
55  int ee = 0;
56  //loop over edges
57  for(;ee<6;ee++) {
58  int jj = 0;
59  if(xH!=NULL) {
60  __CLPK_doublecomplex edge_relative[3];
61  bzero(edge_relative,3*sizeof(__CLPK_doublecomplex));
62  for(;jj<3;jj++) {
63  edge_relative[jj].r = edge_coords[6*ee+jj];
64  edge_relative[jj].r -= edge_coords[6*ee+3+jj];
65  }
66  __CLPK_doublecomplex tmp1 = {1.,0.},tmp2 = {0.,0.};
67  cblas_zgemv(CblasRowMajor,CblasNoTrans,3,3,&tmp1,xH,3,edge_relative,1,&tmp2,&dofs_egdes_X[3*ee],1);
68  }
69  else {
70  for(;jj<3;jj++) {
71  dofs_egdes_X[3*ee+jj].r = edge_coords[6*ee+jj];
72  dofs_egdes_X[3*ee+jj].r -= edge_coords[6*ee+3+jj];
73  }
74  }
75  jj = 0;
76  for(;jj<3;jj++) {
77  lrms += cpow(dofs_egdes_X[3*ee+jj].r+I*dofs_egdes_X[3*ee+jj].i,2.);
78  }
79  }
80  lrms = csqrt( (1./6.)*lrms );
81  (*xlrms).r = creal(lrms);
82  (*xlrms).i = cimag(lrms);
83  PetscFunctionReturn(0);
84 }
void cblas_zgemv(const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA, const int M, const int N, const void *alpha, const void *A, const int lda, const void *X, const int incX, const void *beta, void *Y, const int incY)
Definition: cblas_zgemv.c:12
__CLPK_doublereal r
Definition: lapack_wrap.h:45

◆ calculate_xQ()

PetscErrorCode calculate_xQ ( __CLPK_doublecomplex xlrms,
__CLPK_doublecomplex inv_xH,
__CLPK_doublecomplex dofs_egdes_X_1,
__CLPK_doublecomplex dofs_egdes_X_2,
__CLPK_doublecomplex xQ 
)

Definition at line 85 of file quality_tetrahedorn.c.

87  {
88  PetscFunctionBegin;
89  //add transpose inv_xH
90  int nn = 0;
91  for(;nn<3;nn++) {
92  int mm = 0;
93  for(;mm<3;mm++) {
94  xQ[nn*3+mm].r = inv_xH[mm*3+nn].r;
95  xQ[nn*3+mm].i = inv_xH[mm*3+nn].i;
96  }
97  }
98  //loop over edges
99  double complex tmp0 = -1./(2.*cpow(((*xlrms).r+I*(*xlrms).i),2.));
100  int ee = 0;
101  for(;ee<6;ee++) {
102  int nn = 0;
103  for(;nn<3;nn++) {
104  complex double tmp1 = dofs_egdes_X_1[3*ee+nn].r+I*dofs_egdes_X_1[3*ee+nn].i;
105  int mm = 0;
106  for(;mm<3;mm++) {
107  complex double tmp2 = tmp0*tmp1*(dofs_egdes_X_2[3*ee+mm].r+I*dofs_egdes_X_2[3*ee+mm].i);
108  xQ[mm*3+nn].r += creal(tmp2);
109  xQ[mm*3+nn].i += cimag(tmp2);
110  }
111  }
112  }
113  PetscFunctionReturn(0);
114 }
__CLPK_doublereal i
Definition: lapack_wrap.h:45
__CLPK_doublereal r
Definition: lapack_wrap.h:45

◆ get_edges_from_elem_coords()

PetscErrorCode get_edges_from_elem_coords ( double coords,
double coords_edges 
)

Definition at line 17 of file quality_tetrahedorn.c.

17  {
18  PetscFunctionBegin;
19  cblas_dcopy(3,&coords[0*3],1,&coords_edges[0* 3*2+0],1);
20  cblas_dcopy(3,&coords[1*3],1,&coords_edges[0* 3*2+3],1);
21  cblas_dcopy(3,&coords[0*3],1,&coords_edges[1* 3*2+0],1);
22  cblas_dcopy(3,&coords[2*3],1,&coords_edges[1* 3*2+3],1);
23  cblas_dcopy(3,&coords[0*3],1,&coords_edges[2* 3*2+0],1);
24  cblas_dcopy(3,&coords[3*3],1,&coords_edges[2* 3*2+3],1);
25  cblas_dcopy(3,&coords[1*3],1,&coords_edges[3* 3*2+0],1);
26  cblas_dcopy(3,&coords[2*3],1,&coords_edges[3* 3*2+3],1);
27  cblas_dcopy(3,&coords[1*3],1,&coords_edges[4* 3*2+0],1);
28  cblas_dcopy(3,&coords[3*3],1,&coords_edges[4* 3*2+3],1);
29  cblas_dcopy(3,&coords[2*3],1,&coords_edges[5* 3*2+0],1);
30  cblas_dcopy(3,&coords[3*3],1,&coords_edges[5* 3*2+3],1);
31  PetscFunctionReturn(0);
32 }
void cblas_dcopy(const int N, const double *X, const int incX, double *Y, const int incY)
Definition: cblas_dcopy.c:11

◆ get_qual_ver()

int get_qual_ver ( )

Definition at line 117 of file quality_tetrahedorn.c.

117 { return qual_ver; }
static int qual_ver

◆ quality_volume_length_F()

PetscErrorCode quality_volume_length_F ( double  V,
double alpha2,
double  gamma,
double diffN,
double coords_edges,
double dofs_X,
double dofs_x,
double dofs_iX,
double dofs_ix,
double quality0,
double quality,
double b,
double F,
double iF 
)

Definition at line 154 of file quality_tetrahedorn.c.

156  {
157  PetscFunctionBegin;
158  if(F!=NULL) bzero(F,sizeof(double)*12);
159  if(iF!=NULL) bzero(iF,sizeof(double)*12);
160  double N4[4*4];
162  double H[9];
163  ierr = GradientOfDeformation(diffN,dofs_X,H); CHKERRQ(ierr);
164  double iH[9];
165  if(dofs_iX != NULL) {
166  ierr = GradientOfDeformation(diffN,dofs_iX,iH); CHKERRQ(ierr);
167  } else {
168  bzero(iH,9*sizeof(double));
169  }
170  __CLPK_doublecomplex xH[9];
171  ierr = MakeComplexTensor(H,iH,xH); CHKERRQ(ierr);
172  __CLPK_doublecomplex det_xH;
174  ierr = MakeComplexTensor(H,iH,xH); CHKERRQ(ierr);
175  __CLPK_doublecomplex inv_xH[9];
176  cblas_zcopy(9,xH,1,inv_xH,1);
178  __CLPK_doublecomplex xlrms0;
179  __CLPK_doublecomplex dofs_egdes_Chi[3*6];
180  ierr = calculate_push_edge_relative(coords_edges,NULL,dofs_egdes_Chi,&xlrms0); CHKERRQ(ierr);
182  //printf("%12.10e %12.10e %12.10e\n",xlrms.r-xlrms0.r,xb.r-1,det_xH.r-1);
183  //print_mat_complex(xQ,3,3);
184  //print_mat_complex(xH,3,3);
185  //print_mat_complex(inv_xH,3,3);
186  *quality0 = 6.*sqrt(2.)*(V)/pow(xlrms0.r,3);
187  *quality = xq.r;
188  *b = xb.r;
189  if( F==NULL ) {
190  PetscFunctionReturn(0);
191  }
192  double reQ[9];
193  TakeRe(xQ,reQ);
194  double imQ[9];
195  TakeIm(xQ,imQ);
196  int gg = 0;
197  for(;gg<4;gg++) {
198  double alpha2_val = cblas_ddot(4,&N4[4*gg],1,alpha2,1);
199  int node = 0;
200  for(;node<4;node++) {
201  if(F!=NULL) {
202  F[3*node + 0] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&reQ[0],1);
203  F[3*node + 1] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&reQ[3],1);
204  F[3*node + 2] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&reQ[6],1);
205  }
206  if(iF!=NULL) {
207  iF[3*node + 0] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&imQ[0],1);
208  iF[3*node + 1] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&imQ[3],1);
209  iF[3*node + 2] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&imQ[6],1);
210  }
211  }
212  }
213  PetscFunctionReturn(0);
214 }
static const double G_TET_X4[]
Definition: fem_tools.h:498
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:291
PetscErrorCode TakeRe(__CLPK_doublecomplex *xA, double *reA)
#define QUALITY_VOLUME_LENGTH
PetscErrorCode MakeComplexTensor(double *reA, double *imA, __CLPK_doublecomplex *xA)
Definition: fem_tools.c:428
static PetscErrorCode ierr
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
static const double G_TET_Z4[]
Definition: fem_tools.h:506
PetscErrorCode GradientOfDeformation(double *diffN, double *dofs, double *F)
calculate gradient of deformation
Definition: fem_tools.c:380
CHKERRQ(ierr)
PetscErrorCode TakeIm(__CLPK_doublecomplex *xA, double *imA)
PetscErrorCode DeterminantComplexGradient(__CLPK_doublecomplex *xF, __CLPK_doublecomplex *det_xF)
Definition: fem_tools.c:462
void cblas_zcopy(const int N, const void *X, const int incX, void *Y, const int incY)
Definition: cblas_zcopy.c:11
__CLPK_doublereal r
Definition: lapack_wrap.h:45
static const double G_TET_W4[]
Definition: fem_tools.h:510
PetscErrorCode InvertComplexGradient(__CLPK_doublecomplex *xF)
Definition: fem_tools.c:439
static const double G_TET_Y4[]
Definition: fem_tools.h:502
PetscErrorCode calculate_push_edge_relative(double *edge_coords, __CLPK_doublecomplex *xH, __CLPK_doublecomplex *dofs_egdes_X, __CLPK_doublecomplex *xlrms)

◆ quality_volume_length_K()

int quality_volume_length_K ( double  eps,
double  V,
double alpha2,
double  gamma,
double diffN,
double coords_edges,
double dofs_X,
double dofs_x,
double K,
double Koff 
)

Definition at line 215 of file quality_tetrahedorn.c.

215  {
216  double N4[4*4];
218  double H[9];
219  ierr = GradientOfDeformation(diffN,dofs_X,H); CHKERRQ(ierr);
220  double ZERO[9];
221  bzero(ZERO,sizeof(double)*9);
222  __CLPK_doublecomplex xlrms0;
223  __CLPK_doublecomplex dofs_egdes_Chi[3*6];
224  bzero(K,sizeof(double)*12*12);
225  ierr = calculate_push_edge_relative(coords_edges,NULL,dofs_egdes_Chi,&xlrms0); CHKERRQ(ierr);
226  double _idofs_X[12],_iH[9];
227  int dd = 0;
228  for(;dd<12;dd++) {
229  bzero(_idofs_X,sizeof(double)*12);
230  _idofs_X[dd] = eps;
231  ierr = GradientOfDeformation(diffN,_idofs_X,_iH); CHKERRQ(ierr);
232  __CLPK_doublecomplex xH[9];
233  ierr = MakeComplexTensor(H,_iH,xH); CHKERRQ(ierr);
234  __CLPK_doublecomplex det_xH;
236  ierr = MakeComplexTensor(H,_iH,xH); CHKERRQ(ierr);
237  __CLPK_doublecomplex inv_xH[9];
238  cblas_zcopy(9,xH,1,inv_xH,1);
241  double imQ[9];
242  TakeIm(xQ,imQ);
243  cblas_dscal(9,1./eps,imQ,1);
244  int gg = 0;
245  for(;gg<4;gg++) {
246  double alpha2_val = cblas_ddot(4,&N4[4*gg],1,alpha2,1);
247  int node = 0;
248  for(;node<4;node++) {
249  K[3*12*node + 0*12 + dd] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&imQ[0],1);
250  K[3*12*node + 1*12 + dd] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&imQ[3],1);
251  K[3*12*node + 2*12 + dd] += G_TET_W4[gg]*alpha2_val*cblas_ddot(3,&diffN[node*3+0],1,&imQ[6],1);
252  }
253  }
254  }
255  return 0;
256 }
static const double G_TET_X4[]
Definition: fem_tools.h:498
PetscErrorCode ShapeMBTET(double *N, const double *G_X, const double *G_Y, const double *G_Z, int DIM)
calculate shape functions
Definition: fem_tools.c:291
#define QUALITY_VOLUME_LENGTH
PetscErrorCode MakeComplexTensor(double *reA, double *imA, __CLPK_doublecomplex *xA)
Definition: fem_tools.c:428
static PetscErrorCode ierr
void cblas_dscal(const int N, const double alpha, double *X, const int incX)
Definition: cblas_dscal.c:11
double cblas_ddot(const int N, const double *X, const int incX, const double *Y, const int incY)
Definition: cblas_ddot.c:12
static const double G_TET_Z4[]
Definition: fem_tools.h:506
PetscErrorCode GradientOfDeformation(double *diffN, double *dofs, double *F)
calculate gradient of deformation
Definition: fem_tools.c:380
CHKERRQ(ierr)
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T *> &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:28
PetscErrorCode TakeIm(__CLPK_doublecomplex *xA, double *imA)
PetscErrorCode DeterminantComplexGradient(__CLPK_doublecomplex *xF, __CLPK_doublecomplex *det_xF)
Definition: fem_tools.c:462
void cblas_zcopy(const int N, const void *X, const int incX, void *Y, const int incY)
Definition: cblas_zcopy.c:11
static const double G_TET_W4[]
Definition: fem_tools.h:510
static const double eps
PetscErrorCode InvertComplexGradient(__CLPK_doublecomplex *xF)
Definition: fem_tools.c:439
static const double G_TET_Y4[]
Definition: fem_tools.h:502
PetscErrorCode calculate_push_edge_relative(double *edge_coords, __CLPK_doublecomplex *xH, __CLPK_doublecomplex *dofs_egdes_X, __CLPK_doublecomplex *xlrms)

◆ set_qual_ver()

void set_qual_ver ( int  ver)

Definition at line 116 of file quality_tetrahedorn.c.

116 { qual_ver = ver; }
static int qual_ver

Variable Documentation

◆ ierr

PetscErrorCode ierr
static

Definition at line 15 of file quality_tetrahedorn.c.

◆ qual_ver

int qual_ver = 1
static

Definition at line 115 of file quality_tetrahedorn.c.