source: issm/trunk/src/c/objects/Gauss/GaussTria.cpp@ 12706

/*!\file GaussTria.c
 * \brief: implementation of the GaussTria object
 */

/*Include files: {{{*/
#include "./../objects.h"
/*}}}*/

/*GaussTria constructors and destructors:*/
/*FUNCTION GaussTria::GaussTria() {{{*/
GaussTria::GaussTria(){

        numgauss=-1;

        weights=NULL;
        coords1=NULL;
        coords2=NULL;
        coords3=NULL;

        weight=UNDEF;
        coord1=UNDEF;
        coord2=UNDEF;
        coord3=UNDEF;
}
/*}}}*/
/*FUNCTION GaussTria::GaussTria(int order) {{{*/
GaussTria::GaussTria(int order){

        /*Get gauss points*/
        GaussLegendreTria(&numgauss,&coords1,&coords2,&coords3,&weights,order);

        /*Initialize static fields as undefinite*/
        weight=UNDEF;
        coord1=UNDEF;
        coord2=UNDEF;
        coord3=UNDEF;

}
/*}}}*/
/*FUNCTION GaussTria::GaussTria(int index1,int index2,int order) {{{*/
GaussTria::GaussTria(int index1,int index2,int order){

        /*Intermediaties*/
        IssmPDouble *seg_coords  = NULL;
        IssmPDouble *seg_weights = NULL;
        int     i,index3;

        /*Get Segment gauss points*/
        numgauss=order;
        GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);

        /*Allocate GaussTria fields*/
        coords1=xNew<IssmPDouble>(numgauss);
        coords2=xNew<IssmPDouble>(numgauss);
        coords3=xNew<IssmPDouble>(numgauss);
        weights=xNew<IssmPDouble>(numgauss);

        /*Reverse index1 and 2 if necessary*/
        if (index1>index2){
                index3=index1; index1=index2; index2=index3;
                for(i=0;i<numgauss;i++) seg_coords[i]=-seg_coords[i];
        }

        /*Build Triangle Gauss point*/
        if (index1==0 && index2==1){
                for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
                for(i=0;i<numgauss;i++) coords2[i]=1-0.5*(1.-seg_coords[i]);
                for(i=0;i<numgauss;i++) coords3[i]=0;
                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
        }
        else if (index1==0 && index2==2){
                for(i=0;i<numgauss;i++) coords1[i]=  0.5*(1-seg_coords[i]);
                for(i=0;i<numgauss;i++) coords2[i]= 0 ;
                for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
        }
        else if (index1==1 && index2==2){
                for(i=0;i<numgauss;i++) coords1[i]=0;
                for(i=0;i<numgauss;i++) coords2[i]=  0.5*(1-seg_coords[i]);
                for(i=0;i<numgauss;i++) coords3[i]=1-0.5*(1.-seg_coords[i]);
                for(i=0;i<numgauss;i++) weights[i]=seg_weights[i];
        }
        else
         _error2_("The 2 indices provided are not supported yet (user provided " << index1 << " and " << index2 << ")");

        /*Initialize static fields as undefined*/
        weight=UNDEF;
        coord1=UNDEF;
        coord2=UNDEF;
        coord3=UNDEF;

        /*clean up*/
        xDelete<double>(seg_coords);
        xDelete<double>(seg_weights);
}
/*}}}*/
/*FUNCTION GaussTria::~GaussTria(){{{*/
GaussTria::~GaussTria(){
        xDelete<IssmPDouble>(weights);
        xDelete<IssmPDouble>(coords1);
        xDelete<IssmPDouble>(coords2);
        xDelete<IssmPDouble>(coords3);
}
/*}}}*/

/*Methods*/
/*FUNCTION GaussTria::Echo{{{*/
void GaussTria::Echo(void){

        _printLine_("GaussTria:");
        _printLine_("   numgauss: " << numgauss);

        if (weights){
         _printString_("   weights = ["); 
         for(int i=0;i<numgauss;i++) _printLine_(" " << weights[i]);
         _printLine_("]");
        }
        else _printLine_("weights = NULL");
        if (coords1){
         _printString_("   coords1 = ["); 
         for(int i=0;i<numgauss;i++) _printLine_(" " << coords1[i]);
         _printLine_("]");
        }
        else _printLine_("coords1 = NULL");
        if (coords2){
         _printString_("   coords2 = ["); 
         for(int i=0;i<numgauss;i++) _printLine_(" " << coords2[i]);
         _printLine_("]");
        }
        else _printLine_("coords2 = NULL");
        if (coords3){
         _printString_("   coords3 = ["); 
         for(int i=0;i<numgauss;i++) _printLine_(" " << coords3[i]);
         _printLine_("]");
        }
        else _printLine_("coords3 = NULL");

        _printLine_("   weight = " << weight);
        _printLine_("   coord1 = " << coord1);
        _printLine_("   coord2 = " << coord2);
        _printLine_("   coord3 = " << coord3);

}
/*}}}*/
/*FUNCTION GaussTria::GaussCenter{{{*/
void GaussTria::GaussCenter(void){

        /*update static arrays*/
        coord1=ONETHIRD;
        coord2=ONETHIRD;
        coord3=ONETHIRD;

}
/*}}}*/
/*FUNCTION GaussTria::GaussEdgeCenter{{{*/
void GaussTria::GaussEdgeCenter(int index1,int index2){

        int     index3;

        /*Reverse index1 and 2 if necessary*/
        if (index1>index2){
                index3=index1; index1=index2; index2=index3;
        }

        /*update static arrays*/
        if (index1==0 && index2==1){
                coord1=0.5;
                coord2=0.5;
                coord3=0.0;
        }
        else if (index1==0 && index2==2){
                coord1=0.5;
                coord2=0.0;
                coord3=0.5;
        }
        else if (index1==1 && index2==2){
                coord1=0.0;
                coord2=0.5;
                coord3=0.5;
        }
        else
         _error2_("The 2 indices provided are not supported yet (user provided " << index1 << " and " << index2 << ")");

}
/*}}}*/
/*FUNCTION GaussTria::GaussPoint{{{*/
void GaussTria::GaussPoint(int ig){

        /*Check input in debugging mode*/
         _assert_(ig>=0 && ig< numgauss);

         /*update static arrays*/
         weight=weights[ig];
         coord1=coords1[ig];
         coord2=coords2[ig];
         coord3=coords3[ig];

}
/*}}}*/
/*FUNCTION GaussTria::GaussFromCoords{{{*/
void GaussTria::GaussFromCoords(IssmPDouble x,IssmPDouble y,IssmPDouble* xyz_list){

        /*Intermediaries*/
        IssmPDouble    area = 0;
        IssmPDouble    x1,y1,x2,y2,x3,y3;

        /*in debugging mode: check that the default constructor has been called*/
        _assert_(numgauss==-1);

        x1=*(xyz_list+3*0+0); y1=*(xyz_list+3*0+1);
        x2=*(xyz_list+3*1+0); y2=*(xyz_list+3*1+1);
        x3=*(xyz_list+3*2+0); y3=*(xyz_list+3*2+1);

        area=(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;

        /*Get first area coordinate = det(x-x3  x2-x3 ; y-y3   y2-y3)/area*/
        coord1=((x-x3)*(y2-y3)-(x2-x3)*(y-y3))/area;

        /*Get second area coordinate = det(x1-x3  x-x3 ; y1-y3   y-y3)/area*/
        coord2=((x1-x3)*(y-y3)-(x-x3)*(y1-y3))/area;

        /*Get third  area coordinate 1-area1-area2: */
        coord3=1-coord1-coord2;

}
/*}}}*/
/*FUNCTION GaussTria::GaussVertex{{{*/
void GaussTria::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;
                        break;
                case 1:
                        coord1=0; coord2=1; coord3=0;
                        break;
                case 2:
                        coord1=0; coord2=0; coord3=1;
                        break;
                default:
                        _error2_("vertex index should be in [0 2]");

        }

}
/*}}}*/
/*FUNCTION GaussTria::begin{{{*/
int GaussTria::begin(void){

        /*Check that this has been initialized*/
        _assert_(numgauss>0);
        _assert_(weights);
        _assert_(coords1);
        _assert_(coords2);
        _assert_(coords3);

        /*return first gauss index*/
        return 0;
}
/*}}}*/
/*FUNCTION GaussTria::end{{{*/
int GaussTria::end(void){

        /*Check that this has been initialized*/
        _assert_(numgauss>0);
        _assert_(weights);
        _assert_(coords1);
        _assert_(coords2);
        _assert_(coords3);

        /*return last gauss index +1*/
        return numgauss;
}
/*}}}*/