/*!\file GaussPenta.c * \brief: implementation of the GaussPenta object */ /*Include files: {{{1*/ #include "./../objects.h" /*}}}*/ /*GaussPenta constructors and destructors:*/ /*FUNCTION GaussPenta::GaussPenta() {{{1*/ GaussPenta::GaussPenta(){ numgauss=-1; weights=NULL; coords1=NULL; coords2=NULL; coords3=NULL; coords4=NULL; weight=UNDEF; coord1=UNDEF; coord2=UNDEF; coord3=UNDEF; coord4=UNDEF; } /*}}}*/ /*FUNCTION GaussPenta::GaussPenta(int order_horiz,int order_vert) {{{1*/ GaussPenta::GaussPenta(int order_horiz,int order_vert){ /*Intermediaries*/ int ighoriz,igvert; int numgauss_horiz; int numgauss_vert; double *coords1_horiz = NULL; double *coords2_horiz = NULL; double *coords3_horiz = NULL; double *weights_horiz = NULL; double *coords_vert = NULL; double *weights_vert = NULL; /*Get gauss points*/ GaussLegendreTria(&numgauss_horiz,&coords1_horiz,&coords2_horiz,&coords3_horiz,&weights_horiz,order_horiz); GaussLegendreLinear(&coords_vert,&weights_vert,order_vert); numgauss_vert=order_vert; /*Allocate GaussPenta fields*/ numgauss=numgauss_horiz*numgauss_vert; coords1=(double*)xmalloc(numgauss*sizeof(double)); coords2=(double*)xmalloc(numgauss*sizeof(double)); coords3=(double*)xmalloc(numgauss*sizeof(double)); coords4=(double*)xmalloc(numgauss*sizeof(double)); weights=(double*)xmalloc(numgauss*sizeof(double)); /*Combine Horizontal and vertical points*/ for (ighoriz=0; ighoriz=0 && ig< numgauss); /*update static arrays*/ weight=weights[ig]; coord1=coords1[ig]; coord2=coords2[ig]; coord3=coords3[ig]; coord4=coords4[ig]; } /*}}}*/ /*FUNCTION GaussPenta::GaussVertex{{{1*/ void GaussPenta::GaussVertex(int iv){ /*in debugging mode: check that the default constructor has been called*/ _assert_(numgauss==-1); /*update static arrays*/ switch(iv){ case 0: coord1=1; coord2=0; coord3=0; coord4= -1; break; case 1: coord1=0; coord2=1; coord3=0; coord4= -1; break; case 2: coord1=0; coord2=0; coord3=1; coord4= -1; break; case 3: coord1=1; coord2=0; coord3=0; coord4= +1; break; case 4: coord1=0; coord2=1; coord3=0; coord4= +1; break; case 5: coord1=0; coord2=0; coord3=1; coord4= +1; break; default: _error_("vertex index should be in [0 5]"); } } /*}}}*/ /*FUNCTION GaussPenta::GaussFaceTria{{{1*/ void GaussPenta::GaussFaceTria(int index1, int index2, int index3, int order){ /*in debugging mode: check that the default constructor has been called*/ _assert_(numgauss==-1); /*Basal Tria*/ if(index1==0 && index2==1 && index3==2){ GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order); coords4=(double*)xmalloc(numgauss*sizeof(double)); for(int i=0;i0); _assert_(weights); _assert_(coords1); _assert_(coords2); _assert_(coords3); _assert_(coords4); /*return first gauss index*/ return 0; } /*}}}*/ /*FUNCTION GaussPenta::end{{{1*/ int GaussPenta::end(void){ /*Check that this has been initialized*/ _assert_(numgauss>0); _assert_(weights); _assert_(coords1); _assert_(coords2); _assert_(coords3); _assert_(coords4); /*return last gauss index +1*/ return numgauss; } /*}}}*/ /*FUNCTION GaussPenta::SynchronizeGaussTria{{{1*/ void GaussPenta::SynchronizeGaussTria(GaussTria* gauss_tria){ gauss_tria->coord1=this->coord1; gauss_tria->coord2=this->coord2; gauss_tria->coord3=this->coord3; gauss_tria->weight=UNDEF; } /*}}}*/