Ignore:
Timestamp:
09/16/13 09:43:55 (12 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 16135

Location:
issm/trunk
Files:
1 deleted
9 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        1414probe.results
        1515stXXXX*
        16 
         16.deps
         17.dirstamp

  • issm/trunk/src/c/classes

    • Property svn:ignore set to
      .deps
      .dirstamp

  • issm/trunk/src/c/classes/gauss

    • Property svn:ignore set to
      .deps
      .dirstamp

  • issm/trunk/src/c/classes/gauss/GaussPenta.cpp

    r15396 r16137  
    88#include "../../shared/Exceptions/exceptions.h"
    99#include "../../shared/MemOps/MemOps.h"
     10#include "../../shared/Enum/Enum.h"
    1011#include "../../shared/Numerics/GaussPoints.h"
    1112#include "../../shared/Numerics/constants.h"
    … …  
    237238}
    238239/*}}}*/
     240/*FUNCTION GaussPenta::GaussPenta(int index,double r1,double r2,int order) {{{*/
     241GaussPenta::GaussPenta(int index,IssmDouble r1,IssmDouble r2,bool mainlyfloating,int order){
     242
     243        /*
     244         *  ^
     245         *  |
     246         * 1|\
     247         *  |  \
     248         *  |    \
     249         *  |      \
     250         *  |        \
     251         *  |          \
     252         *  |    +(x,y)  \
     253         *  |              \
     254         *  +---------------+-->
     255         *  0               1
     256         *
     257         */
     258        int         ig;
     259        IssmDouble x,y;
     260        IssmDouble xy_list[3][2];
     261
     262        if(mainlyfloating){
     263                /*Get gauss points*/
     264                GaussLegendreTria(&this->numgauss,&this->coords1,&this->coords2,&this->coords3,&this->weights,order);
     265
     266                xy_list[0][0]=0;  xy_list[0][1]=0;
     267                xy_list[1][0]=r1; xy_list[1][1]=0;
     268                xy_list[2][0]=0;  xy_list[2][1]=r2;
     269
     270                for(ig=0;ig<this->numgauss;ig++){
     271                        x = this->coords1[ig]*xy_list[0][0] + this->coords2[ig]*xy_list[1][0] + this->coords3[ig]*xy_list[2][0];
     272                        y = this->coords1[ig]*xy_list[0][1] + this->coords2[ig]*xy_list[1][1] + this->coords3[ig]*xy_list[2][1];
     273
     274                        switch(index){
     275                                case 0:
     276                                        this->coords1[ig] = 1.-x-y;
     277                                        this->coords2[ig] = x;
     278                                        this->coords3[ig] = y;
     279                                        break;
     280                                case 1:
     281                                        this->coords1[ig] = y;
     282                                        this->coords2[ig] = 1.-x-y;
     283                                        this->coords3[ig] = x;
     284                                        break;
     285                                case 2:
     286                                        this->coords1[ig] = x;
     287                                        this->coords2[ig] = y;
     288                                        this->coords3[ig] = 1.-x-y;
     289                                        break;
     290                                default:
     291                                        _error_("index "<<index<<" not supported yet");
     292                        }
     293                        this->weights[ig] = this->weights[ig]*r1*r2;
     294                        this->coords4=xNew<IssmDouble>(numgauss);
     295                        for(ig=0;ig<numgauss;ig++) this->coords4[ig]=-1.0;
     296                }
     297        }
     298        else{
     299                /*Double number of gauss points*/
     300                GaussPenta *gauss1    = NULL;
     301                GaussPenta *gauss2    = NULL;
     302                gauss1=new GaussPenta(0,1,2,order);
     303                gauss2=new GaussPenta(0,1,2,order);
     304
     305                xy_list[0][0]=r1; xy_list[0][1]=0;
     306                xy_list[1][0]=0;  xy_list[1][1]=1.;
     307                xy_list[2][0]=0;  xy_list[2][1]=r2;
     308
     309                        //gauss1->Echo();
     310                for(ig=0;ig<gauss1->numgauss;ig++){
     311                        x = gauss1->coords1[ig]*xy_list[0][0] + gauss1->coords2[ig]*xy_list[1][0] + gauss1->coords3[ig]*xy_list[2][0];
     312                        y = gauss1->coords1[ig]*xy_list[0][1] + gauss1->coords2[ig]*xy_list[1][1] + gauss1->coords3[ig]*xy_list[2][1];
     313
     314                        switch(index){
     315                                case 0:
     316                                        gauss1->coords1[ig] = 1.-x-y;
     317                                        gauss1->coords2[ig] = x;
     318                                        gauss1->coords3[ig] = y;
     319                                        break;
     320                                case 1:
     321                                        gauss1->coords1[ig] = y;
     322                                        gauss1->coords2[ig] = 1.-x-y;
     323                                        gauss1->coords3[ig] = x;
     324                                        break;
     325                                case 2:
     326                                        gauss1->coords1[ig] = x;
     327                                        gauss1->coords2[ig] = y;
     328                                        gauss1->coords3[ig] = 1.-x-y;
     329                                        break;
     330                                default:
     331                                        _error_("index "<<index<<" not supported yet");
     332                        }
     333                        gauss1->weights[ig] = gauss1->weights[ig]*r1*(1-r2);
     334                }
     335                        //gauss1->Echo();
     336                xy_list[0][0]=r1; xy_list[0][1]=0;
     337                xy_list[1][0]=1.; xy_list[1][1]=0;
     338                xy_list[2][0]=0;  xy_list[2][1]=1.;
     339
     340                        //gauss2->Echo();
     341                for(ig=0;ig<gauss2->numgauss;ig++){
     342                        x = gauss2->coords1[ig]*xy_list[0][0] + gauss2->coords2[ig]*xy_list[1][0] + gauss2->coords3[ig]*xy_list[2][0];
     343                        y = gauss2->coords1[ig]*xy_list[0][1] + gauss2->coords2[ig]*xy_list[1][1] + gauss2->coords3[ig]*xy_list[2][1];
     344
     345                        switch(index){
     346                                case 0:
     347                                        gauss2->coords1[ig] = 1.-x-y;
     348                                        gauss2->coords2[ig] = x;
     349                                        gauss2->coords3[ig] = y;
     350                                        break;
     351                                case 1:
     352                                        gauss2->coords1[ig] = y;
     353                                        gauss2->coords2[ig] = 1.-x-y;
     354                                        gauss2->coords3[ig] = x;
     355                                        break;
     356                                case 2:
     357                                        gauss2->coords1[ig] = x;
     358                                        gauss2->coords2[ig] = y;
     359                                        gauss2->coords3[ig] = 1.-x-y;
     360                                        break;
     361                                default:
     362                                        _error_("index "<<index<<" not supported yet");
     363                        }
     364                        gauss2->weights[ig] = gauss2->weights[ig]*(1-r1);
     365                }
     366
     367                this->numgauss = gauss1->numgauss + gauss2->numgauss;
     368                this->coords1=xNew<IssmDouble>(this->numgauss);
     369                this->coords2=xNew<IssmDouble>(this->numgauss);
     370                this->coords3=xNew<IssmDouble>(this->numgauss);
     371                this->coords4=xNew<IssmDouble>(this->numgauss);
     372                this->weights=xNew<IssmDouble>(this->numgauss);
     373
     374                for(ig=0;ig<gauss1->numgauss;ig++){ // Add the first triangle gauss points
     375                        this->coords1[ig]=gauss1->coords1[ig];
     376                        this->coords2[ig]=gauss1->coords2[ig];
     377                        this->coords3[ig]=gauss1->coords3[ig];
     378                        this->coords4[ig]=gauss1->coords4[ig];
     379                        this->weights[ig]=gauss1->weights[ig];
     380                }
     381                for(ig=0;ig<gauss2->numgauss;ig++){ // Add the second triangle gauss points
     382                        this->coords1[gauss1->numgauss+ig]=gauss2->coords1[ig];
     383                        this->coords2[gauss1->numgauss+ig]=gauss2->coords2[ig];
     384                        this->coords3[gauss1->numgauss+ig]=gauss2->coords3[ig];
     385                        this->coords4[gauss1->numgauss+ig]=gauss2->coords4[ig];
     386                        this->weights[gauss1->numgauss+ig]=gauss2->weights[ig];
     387                }
     388
     389                /*Delete gauss points*/
     390                delete gauss1;
     391                delete gauss2;
     392        }
     393
     394        /*Initialize static fields as undefined*/
     395        weight=UNDEF;
     396        coord1=UNDEF;
     397        coord2=UNDEF;
     398        coord3=UNDEF;
     399        coord4=UNDEF;
     400}
     401/*}}}*/
     402/*FUNCTION GaussPenta::GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert){{{*/
     403GaussPenta::GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert){
     404
     405        /*Intermediaties*/
     406        IssmPDouble *seg_horiz_coords  = NULL;
     407        IssmPDouble *seg_horiz_weights = NULL;
     408        IssmPDouble *seg_vert_coords   = NULL;
     409        IssmPDouble *seg_vert_weights  = NULL;
     410
     411        /*get the gauss points using the product of two line rules*/
     412        GaussLegendreLinear(&seg_horiz_coords,&seg_horiz_weights,order_horiz);
     413        GaussLegendreLinear(&seg_vert_coords, &seg_vert_weights, order_vert);
     414
     415        /*Allocate GaussPenta fields*/
     416        numgauss=order_horiz*order_vert;
     417        coords1=xNew<IssmDouble>(numgauss);
     418        coords2=xNew<IssmDouble>(numgauss);
     419        coords3=xNew<IssmDouble>(numgauss);
     420        coords4=xNew<IssmDouble>(numgauss);
     421        weights=xNew<IssmDouble>(numgauss);
     422
     423        /*Quads: get the gauss points using the product of two line rules  */
     424        for(int i=0;i<order_horiz;i++){
     425                for(int j=0;j<order_vert;j++){
     426                        coords1[i*order_vert+j]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
     427                        coords2[i*order_vert+j]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
     428                        coords3[i*order_vert+j]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_horiz_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
     429                        coords4[i*order_vert+j]=seg_vert_coords[j];
     430                        weights[i*order_vert+j]=seg_horiz_weights[i]*seg_vert_weights[j];
     431                }
     432        }
     433
     434        /*clean-up*/
     435        xDelete<IssmPDouble>(seg_horiz_coords);
     436        xDelete<IssmPDouble>(seg_horiz_weights);
     437        xDelete<IssmPDouble>(seg_vert_coords);
     438        xDelete<IssmPDouble>(seg_vert_weights);
     439}
     440/*}}}*/
    239441/*FUNCTION GaussPenta::~GaussPenta(){{{*/
    240442GaussPenta::~GaussPenta(){
    … …  
    327529        /*update static arrays*/
    328530        switch(iv){
    329                 case 0:
    330                         coord1=1; coord2=0; coord3=0; coord4= -1;
    331                         break;
    332                 case 1:
    333                         coord1=0; coord2=1; coord3=0; coord4= -1;
    334                         break;
    335                 case 2:
    336                         coord1=0; coord2=0; coord3=1; coord4= -1;
    337                         break;
    338                 case 3:
    339                         coord1=1; coord2=0; coord3=0; coord4= +1;
    340                         break;
    341                 case 4:
    342                         coord1=0; coord2=1; coord3=0; coord4= +1;
    343                         break;
    344                 case 5:
    345                         coord1=0; coord2=0; coord3=1; coord4= +1;
    346                         break;
    347                 default:
    348                         _error_("vertex index should be in [0 5]");
     531                case 0: coord1=1.; coord2=0.; coord3=0.; coord4= -1.; break;
     532                case 1: coord1=0.; coord2=1.; coord3=0.; coord4= -1.; break;
     533                case 2: coord1=0.; coord2=0.; coord3=1.; coord4= -1.; break;
     534                case 3: coord1=1.; coord2=0.; coord3=0.; coord4= +1.; break;
     535                case 4: coord1=0.; coord2=1.; coord3=0.; coord4= +1.; break;
     536                case 5: coord1=0.; coord2=0.; coord3=1.; coord4= +1.; break;
     537                default: _error_("vertex index should be in [0 5]");
    349538
    350539        }
    … …  
    366555        else{
    367556                _error_("Tria not supported yet");
     557        }
     558
     559}
     560/*}}}*/
     561/*FUNCTION GaussPenta::GaussNode{{{*/
     562void GaussPenta::GaussNode(int finiteelement,int iv){
     563
     564        /*in debugging mode: check that the default constructor has been called*/
     565        _assert_(numgauss==-1);
     566
     567        /*update static arrays*/
     568        switch(finiteelement){
     569                case P1Enum: case P1DGEnum:
     570                        switch(iv){
     571                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     572                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     573                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     574                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     575                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     576                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     577                                default: _error_("node index should be in [0 5]");
     578                        }
     579                        break;
     580                case P1xP2Enum:
     581                        switch(iv){
     582                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     583                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     584                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     585                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     586                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     587                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     588
     589                                case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
     590                                case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
     591                                case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
     592                                default: _error_("node index should be in [0 8]");
     593                        }
     594                        break;
     595                case P2xP1Enum:
     596                        switch(iv){
     597                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     598                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     599                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     600                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     601                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     602                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     603
     604                                case  6: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
     605                                case  7: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
     606                                case  8: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
     607                                case  9: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
     608                                case 10: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
     609                                case 11: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
     610                                default: _error_("node index should be in [0 11]");
     611                        }
     612                        break;
     613                case P1bubbleEnum:  case P1bubblecondensedEnum:
     614                        switch(iv){
     615                                case 0: coord1=1.;    coord2=0.;    coord3=0.;    coord4=-1.; break;
     616                                case 1: coord1=0.;    coord2=1.;    coord3=0.;    coord4=-1.; break;
     617                                case 2: coord1=0.;    coord2=0.;    coord3=1.;    coord4=-1.; break;
     618                                case 3: coord1=1.;    coord2=0.;    coord3=0.;    coord4=+1.; break;
     619                                case 4: coord1=0.;    coord2=1.;    coord3=0.;    coord4=+1.; break;
     620                                case 5: coord1=0.;    coord2=0.;    coord3=1.;    coord4=+1.; break;
     621                                case 6: coord1=1./3.; coord2=1./3.; coord3=1./3.; coord4=0.;  break;
     622                                default: _error_("node index should be in [0 6]");
     623                        }
     624                        break;
     625                case P2Enum:
     626                        switch(iv){
     627                                case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
     628                                case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
     629                                case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
     630                                case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
     631                                case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
     632                                case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
     633
     634                                case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
     635                                case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
     636                                case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
     637
     638                                case  9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
     639                                case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
     640                                case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
     641                                case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
     642                                case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
     643                                case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
     644                                default: _error_("node index should be in [0 14]");
     645                        }
     646                        break;
     647                default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
    368648        }
    369649

  • issm/trunk/src/c/classes/gauss/GaussPenta.h

    r15396 r16137  
    3535                GaussPenta(int index1, int index2, int index3, int order);
    3636                GaussPenta(int index1, int index2, int index3, int index4,int order_horiz,int order_vert);
     37                GaussPenta(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
     38                GaussPenta(IssmDouble area_coordinates[4][3],int order_horiz,int order_vert);
    3739                ~GaussPenta();
    3840
    … …  
    4345                void GaussPoint(int ig);
    4446                void GaussVertex(int iv);
     47                void GaussNode(int finitelement,int iv);
    4548                void GaussFaceTria(int index1, int index2, int index3, int order);
    4649                void GaussCenter(void);

  • issm/trunk/src/c/classes/gauss/GaussTria.cpp

    r15396 r16137  
    9999        xDelete<double>(seg_coords);
    100100        xDelete<double>(seg_weights);
     101}
     102/*}}}*/
     103/*FUNCTION GaussTria::GaussTria(IssmDouble area_coordinates,int order) {{{*/
     104GaussTria::GaussTria(IssmDouble area_coordinates[2][3],int order){
     105
     106        /*Intermediaties*/
     107        IssmPDouble *seg_coords  = NULL;
     108        IssmPDouble *seg_weights = NULL;
     109        int     i,index3;
     110
     111        /*Get Segment gauss points*/
     112        numgauss=order;
     113        GaussLegendreLinear(&seg_coords,&seg_weights,numgauss);
     114
     115        /*Allocate GaussTria fields*/
     116        coords1=xNew<IssmDouble>(numgauss);
     117        coords2=xNew<IssmDouble>(numgauss);
     118        coords3=xNew<IssmDouble>(numgauss);
     119        weights=xNew<IssmDouble>(numgauss);
     120
     121        /*Build Triangle Gauss point*/
     122        for(i=0;i<numgauss;i++){
     123                coords1[i]=0.5*(area_coordinates[0][0]+area_coordinates[1][0]) + 0.5*seg_coords[i]*(area_coordinates[1][0]-area_coordinates[0][0]);
     124                coords2[i]=0.5*(area_coordinates[0][1]+area_coordinates[1][1]) + 0.5*seg_coords[i]*(area_coordinates[1][1]-area_coordinates[0][1]);
     125                coords3[i]=0.5*(area_coordinates[0][2]+area_coordinates[1][2]) + 0.5*seg_coords[i]*(area_coordinates[1][2]-area_coordinates[0][2]);
     126                weights[i]=seg_weights[i];
     127        }
     128
     129        /*Initialize static fields as undefined*/
     130        weight=UNDEF;
     131        coord1=UNDEF;
     132        coord2=UNDEF;
     133        coord3=UNDEF;
     134
     135        /*clean up*/
     136        xDelete<IssmPDouble>(seg_coords);
     137        xDelete<IssmPDouble>(seg_weights);
    101138}
    102139/*}}}*/
    … …  
    260297GaussTria::~GaussTria(){
    261298        xDelete<IssmDouble>(weights);
     299        xDelete<IssmDouble>(coords3);
     300        xDelete<IssmDouble>(coords2);
    262301        xDelete<IssmDouble>(coords1);
    263         xDelete<IssmDouble>(coords2);
    264         xDelete<IssmDouble>(coords3);
     302
    265303}
    266304/*}}}*/
    … …  
    404442/*}}}*/
    405443/*FUNCTION GaussTria::GaussNode{{{*/
    406 void GaussTria::GaussNode(int numnodes,int iv){
     444void GaussTria::GaussNode(int finiteelement,int iv){
    407445
    408446        /*in debugging mode: check that the default constructor has been called*/
    … …  
    410448
    411449        /*update static arrays*/
    412         switch(numnodes){
    413                 case 3: //P1 Lagrange element
     450        switch(finiteelement){
     451                case P1Enum: case P1DGEnum:
    414452                        switch(iv){
    415453                                case 0: coord1=1.; coord2=0.; coord3=0.; break;
    … …  
    419457                        }
    420458                        break;
    421                 case 6: //P2 Lagrange element
     459                case P1bubbleEnum: case P1bubblecondensedEnum:
     460                        switch(iv){
     461                                case 0: coord1=1.;    coord2=0.;    coord3=0.;    break;
     462                                case 1: coord1=0.;    coord2=1.;    coord3=0.;    break;
     463                                case 2: coord1=0.;    coord2=0.;    coord3=1.;    break;
     464                                case 3: coord1=1./3.; coord2=1./3.; coord3=1./3.; break;
     465                                default: _error_("node index should be in [0 3]");
     466                        }
     467                        break;
     468                case P2Enum:
    422469                        switch(iv){
    423470                                case 0: coord1=1.; coord2=0.; coord3=0.; break;
    … …  
    430477                        }
    431478                        break;
    432                 default: _error_("supported number of nodes are 3 and 6");
     479                default: _error_("Finite element "<<EnumToStringx(finiteelement)<<" not supported");
    433480        }
    434481

  • issm/trunk/src/c/classes/gauss/GaussTria.h

    r15396 r16137  
    3131                GaussTria(int index1,int index2,int order);
    3232                GaussTria(int index,IssmDouble r1, IssmDouble r2,bool maintlyfloating,int order);
     33                GaussTria(IssmDouble area_coordinates[3][3],int order);
    3334                ~GaussTria();
    3435
    … …  
    4041                void GaussPoint(int ig);
    4142                void GaussVertex(int iv);
    42                 void GaussNode(int numnodes,int iv);
     43                void GaussNode(int finitelement,int iv);
    4344                void GaussCenter(void);
    4445                void GaussEdgeCenter(int index1,int index2);
Note: See TracChangeset for help on using the changeset viewer.