Changeset 6029


Ignore:
Timestamp:
09/24/10 13:23:32 (14 years ago)
Author:
Mathieu Morlighem
Message:

Fixed Adjoint Balanced thickness-> The stiffness matrix must be transposed

Location:
issm/trunk/src/c/objects
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk/src/c/objects/Elements/Tria.cpp

    r5989 r6029  
    718718                        Ke=CreateKMatrixPrognostic();
    719719                        break;
    720                 case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum:
     720                case BalancedthicknessAnalysisEnum:
    721721                        Ke=CreateKMatrixBalancedthickness();
     722                        break;
     723                case AdjointBalancedthicknessAnalysisEnum:
     724                        Ke=CreateKMatrixAdjointBalancedthickness();
    722725                        break;
    723726                case BalancedvelocitiesAnalysisEnum:
    … …  
    23842387
    23852388/*Tria specific routines: */
     2389/*FUNCTION Tria::CreateKMatrixAdjointBalancedthickness {{{1*/
     2390ElementMatrix* Tria::CreateKMatrixAdjointBalancedthickness(void){
     2391
     2392        ElementMatrix* Ke=NULL;
     2393
     2394        /*Get Element Matrix of the forward model*/
     2395        switch(GetElementType()){
     2396                case P1Enum:
     2397                        Ke=CreateKMatrixBalancedthickness_CG();
     2398                        break;
     2399                case P1DGEnum:
     2400                        Ke=CreateKMatrixBalancedthickness_DG();
     2401                        break;
     2402                default:
     2403                        ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType()));
     2404        }
     2405
     2406        /*Transpose and return Ke*/
     2407        Ke->Transpose();
     2408        return Ke;
     2409}
     2410/*}}}*/
    23862411/*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/
    23872412ElementMatrix* Tria::CreateKMatrixBalancedthickness(void){

  • TabularUnified issm/trunk/src/c/objects/Elements/Tria.h

    r5933 r6029  
    121121                /*}}}*/
    122122                /*Tria specific routines:{{{1*/
     123                ElementMatrix* CreateKMatrixAdjointBalancedthickness(void);
    123124                ElementMatrix* CreateKMatrixBalancedthickness(void);
    124125                ElementMatrix* CreateKMatrixBalancedthickness_DG(void);

  • TabularUnified issm/trunk/src/c/objects/Loads/Numericalflux.cpp

    r6026 r6029  
    509509
    510510        /*Initialize Element matrix and return if necessary*/
     511        ElementMatrix* Ke = NULL;
    511512        Tria*  tria=(Tria*)element;
    512513        if(tria->IsOnWater()) return NULL;
    513         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    514514
    515515        /*Retrieve all inputs and parameters*/
    … …  
    532532        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    533533        if (UdotN<=0){
    534                 /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    535                 return NULL;
     534                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     535        }
     536        else{
     537                Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    536538        }
    537539
    … …  
    659661
    660662        /*Initialize Element matrix and return if necessary*/
     663        ElementMatrix* Ke = NULL;
    661664        Tria*  tria=(Tria*)element;
    662665        if(tria->IsOnWater()) return NULL;
    663         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    664666
    665667        /*Retrieve all inputs and parameters*/
    … …  
    681683        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    682684        if (UdotN<=0){
    683                 /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
    684                 return NULL;
     685                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
     686        }
     687        else{
     688                Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    685689        }
    686690
    … …  
    731735ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessInternal(void){
    732736
    733         return CreateKMatrixBalancedthicknessInternal();
     737        ElementMatrix* Ke=CreateKMatrixBalancedthicknessInternal();
     738        if (Ke) Ke->Transpose();
     739        return Ke;
    734740}
    735741/*}}}*/
    … …  
    737743ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessBoundary(void){
    738744
    739         return CreateKMatrixBalancedthicknessBoundary();
     745        ElementMatrix* Ke=CreateKMatrixBalancedthicknessBoundary();
     746        if(Ke) Ke->Transpose();
     747        return Ke;
    740748}
    741749/*}}}*/
    … …  
    779787
    780788        /*Initialize Load Vector and return if necessary*/
     789        ElementVector* pe = NULL;
    781790        Tria*  tria=(Tria*)element;
    782791        if(tria->IsOnWater()) return NULL;
    783         ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    784792
    785793        /*Retrieve all inputs and parameters*/
    … …  
    800808        vyaverage_input->GetParameterValue(&mean_vy,gauss);
    801809        delete gauss;
     810
    802811        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    803812        if (UdotN>0){
    804                 /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    805                 return NULL;
     813                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
     814        }
     815        else{
     816                ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    806817        }
    807818
    … …  
    868879
    869880        /*Initialize Load Vector and return if necessary*/
     881        ElementVector* pe = NULL;
    870882        Tria*  tria=(Tria*)element;
    871883        if(tria->IsOnWater()) return NULL;
    872         ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    873884
    874885        /*Retrieve all inputs and parameters*/
    … …  
    876887        Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input);
    877888        Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input);
    878         Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum);
    879 
    880         /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/
    881         if (!thickness_input){
    882                 thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    883         }
     889        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
    884890        GetNormal(&normal[0],xyz_list);
    885891
    … …  
    895901        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    896902        if (UdotN>0){
    897                 /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    898                 return NULL;
     903                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
     904        }
     905        else{
     906                pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    899907        }
    900908
    … …  
    925933ElementVector* Numericalflux::CreatePVectorAdjointBalancedthickness(void){
    926934
    927         int type;
    928         inputs->GetParameterValue(&type,TypeEnum);
    929 
    930         switch(type){
    931                 case InternalEnum:
    932                         return CreatePVectorAdjointBalancedthicknessInternal();
    933                 case BoundaryEnum:
    934                         return CreatePVectorAdjointBalancedthicknessBoundary();
    935                 default:
    936                         ISSMERROR("type not supported yet");
    937         }
    938 }
    939 /*}}}*/
    940 /*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessInternal{{{1*/
    941 ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessInternal(void){
    942 
    943         /*Nothing added to PVector*/
     935        /*No PVector for the Adjoint*/
    944936        return NULL;
    945 
    946 }
    947 /*}}}*/
    948 /*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary{{{1*/
    949 ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary(void){
    950 
    951         /* constants*/
    952         const int numdof=NDOF1*NUMVERTICES_BOUNDARY;
    953 
    954         /* Intermediaries*/
    955         int        i,j,ig,index1,index2;
    956         double     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
    957         double     xyz_list[NUMVERTICES_BOUNDARY][3];
    958         double     normal[2];
    959         double     L[numdof];
    960         GaussTria *gauss;
    961 
    962         /*Initialize Load Vector and return if necessary*/
    963         Tria*  tria=(Tria*)element;
    964         if(tria->IsOnWater()) return NULL;
    965         ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
    966 
    967         /*Retrieve all inputs and parameters*/
    968         GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
    969         Input* vxaverage_input=tria->inputs->GetInput(VxEnum);           ISSMASSERT(vxaverage_input);
    970         Input* vyaverage_input=tria->inputs->GetInput(VyEnum);           ISSMASSERT(vyaverage_input);
    971         Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); ISSMASSERT(thickness_input);
    972         GetNormal(&normal[0],xyz_list);
    973 
    974         /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
    975         index1=tria->GetNodeIndex(nodes[0]);
    976         index2=tria->GetNodeIndex(nodes[1]);
    977 
    978         gauss=new GaussTria();
    979         gauss->GaussEdgeCenter(index1,index2);
    980         vxaverage_input->GetParameterValue(&mean_vx,gauss);
    981         vyaverage_input->GetParameterValue(&mean_vy,gauss);
    982         delete gauss;
    983         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
    984         if (UdotN>0){
    985                 /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
    986                 return NULL;
    987         }
    988 
    989         /* Start  looping on the number of gaussian points: */
    990         gauss=new GaussTria(index1,index2,2);
    991         for(ig=gauss->begin();ig<gauss->end();ig++){
    992 
    993                 gauss->GaussPoint(ig);
    994 
    995                 tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
    996 
    997                 vxaverage_input->GetParameterValue(&vx,gauss);
    998                 vyaverage_input->GetParameterValue(&vy,gauss);
    999                 thickness_input->GetParameterValue(&thickness,gauss);
    1000                 UdotN=vx*normal[0]+vy*normal[1];
    1001                 tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
    1002                 DL= - gauss->weight*Jdet*UdotN*thickness;
    1003 
    1004                 for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
    1005         }
    1006 
    1007         /*Clean up and return*/
    1008         delete gauss;
    1009         return pe;
    1010937}
    1011938/*}}}*/

  • TabularUnified issm/trunk/src/c/objects/Loads/Numericalflux.h

    r6026 r6029  
    9090                ElementVector* CreatePVectorBalancedthicknessBoundary(void);
    9191                ElementVector* CreatePVectorAdjointBalancedthickness(void);
    92                 ElementVector* CreatePVectorAdjointBalancedthicknessInternal(void);
    93                 ElementVector* CreatePVectorAdjointBalancedthicknessBoundary(void);
    9492                /*}}}*/
    9593

  • TabularUnified issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp

    r6027 r6029  
    4343        this->col_sglobaldoflist=NULL;
    4444
     45}
     46/*}}}*/
     47/*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke){{{1*/
     48ElementMatrix::ElementMatrix(ElementMatrix* Ke){
     49
     50        if(!Ke) ISSMERROR("Input Element Matrix is a NULL pointer");
     51        this->Init(Ke);
     52        return;
    4553}
    4654/*}}}*/
    … …  
    319327
    320328        /*Intermediaries*/
    321         double *values_copy=NULL;
    322         int     temp;
     329        ElementMatrix* Ke_copy=new ElementMatrix(this);
     330
     331        /*Update sizes*/
     332        this->nrows=Ke_copy->ncols;
     333        this->ncols=Ke_copy->nrows;
     334
     335        /*Transpose values*/
     336        for (int i=0;i<this->nrows;i++) for(int j=0;j<this->ncols;j++) this->values[i*this->ncols+j]=Ke_copy->values[j*Ke_copy->ncols+i];
    323337
    324338        /*Transpose indices*/
    … …  
    327341        }
    328342
    329         /*Transpose values*/
    330         values_copy=(double*)xmalloc(this->nrows*this->ncols*sizeof(double));
    331         memcpy(values_copy,this->values,this->nrows*this->ncols*sizeof(double));
    332         for (int i=0;i<this->nrows;i++) for(int j=0;j<this->ncols;j++) this->values[j*this->nrows+i]=values_copy[i*this->ncols+j];
    333 
    334         /*Update sizes*/
    335         temp=this->nrows;
    336         this->nrows=this->ncols;
    337         this->ncols=temp;
    338 
    339343        /*Clean up and return*/
    340         xfree((void**)&values_copy);
     344        delete Ke_copy;
     345        return;
    341346}
    342347/*}}}*/

  • TabularUnified issm/trunk/src/c/objects/Numerics/ElementMatrix.h

    r6027 r6029  
    5252                /*ElementMatrix constructors, destructors {{{1*/
    5353                ElementMatrix();
     54                ElementMatrix(ElementMatrix* Ke);
    5455                ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2);
    5556                ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2,ElementMatrix* Ke3);
Note: See TracChangeset for help on using the changeset viewer.