Changeset 6029 for issm/trunk/src/c/objects/Loads/Numericalflux.cpp

Timestamp:

09/24/10 13:23:32 (14 years ago)

Author:

Mathieu Morlighem

Message:

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

File:

• issm/trunk/src/c/objects/Loads/Numericalflux.cpp (modified) (12 diffs)

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

@@ -509 +509 @@
 
         /*Initialize Element matrix and return if necessary*/
+        ElementMatrix* Ke = NULL;
         Tria*  tria=(Tria*)element;
         if(tria->IsOnWater()) return NULL;
-        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
 
         /*Retrieve all inputs and parameters*/
@@ -532 +532 @@
         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
         if (UdotN<=0){
-                /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
-                return NULL;
+                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
+        }
+        else{
+                Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
         }
 
@@ -659 +661 @@
 
         /*Initialize Element matrix and return if necessary*/
+        ElementMatrix* Ke = NULL;
         Tria*  tria=(Tria*)element;
         if(tria->IsOnWater()) return NULL;
-        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
 
         /*Retrieve all inputs and parameters*/
@@ -681 +683 @@
         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
         if (UdotN<=0){
-                /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
-                return NULL;
+                return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/
+        }
+        else{
+                Ke=new ElementMatrix(nodes,NUMVERTICES_BOUNDARY,this->parameters);
         }
 
@@ -731 +735 @@
 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessInternal(void){
 
-        return CreateKMatrixBalancedthicknessInternal();
+        ElementMatrix* Ke=CreateKMatrixBalancedthicknessInternal();
+        if (Ke) Ke->Transpose();
+        return Ke;
 }
 /*}}}*/
@@ -737 +743 @@
 ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessBoundary(void){
 
-        return CreateKMatrixBalancedthicknessBoundary();
+        ElementMatrix* Ke=CreateKMatrixBalancedthicknessBoundary();
+        if(Ke) Ke->Transpose();
+        return Ke;
 }
 /*}}}*/
@@ -779 +787 @@
 
         /*Initialize Load Vector and return if necessary*/
+        ElementVector* pe = NULL;
         Tria*  tria=(Tria*)element;
         if(tria->IsOnWater()) return NULL;
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
 
         /*Retrieve all inputs and parameters*/
@@ -800 +808 @@
         vyaverage_input->GetParameterValue(&mean_vy,gauss);
         delete gauss;
+
         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
         if (UdotN>0){
-                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
-                return NULL;
+                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
+        }
+        else{
+                ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
         }
 
@@ -868 +879 @@
 
         /*Initialize Load Vector and return if necessary*/
+        ElementVector* pe = NULL;
         Tria*  tria=(Tria*)element;
         if(tria->IsOnWater()) return NULL;
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
 
         /*Retrieve all inputs and parameters*/
@@ -876 +887 @@
         Input* vxaverage_input=tria->inputs->GetInput(VxEnum); ISSMASSERT(vxaverage_input); 
         Input* vyaverage_input=tria->inputs->GetInput(VyEnum); ISSMASSERT(vyaverage_input);
-        Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum);
-
-        /*Here, as it is a forcing, we have H=Hobs by default (for control methods)*/
-        if (!thickness_input){
-                thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
-        }
+        Input* thickness_input=tria->inputs->GetInput(ThicknessEnum); ISSMASSERT(thickness_input);
         GetNormal(&normal[0],xyz_list);
 
@@ -895 +901 @@
         UdotN=mean_vx*normal[0]+mean_vy*normal[1];
         if (UdotN>0){
-                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
-                return NULL;
+                return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
+        }
+        else{
+                pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
         }
 
@@ -925 +933 @@
 ElementVector* Numericalflux::CreatePVectorAdjointBalancedthickness(void){
 
-        int type;
-        inputs->GetParameterValue(&type,TypeEnum);
-
-        switch(type){
-                case InternalEnum:
-                        return CreatePVectorAdjointBalancedthicknessInternal();
-                case BoundaryEnum:
-                        return CreatePVectorAdjointBalancedthicknessBoundary();
-                default:
-                        ISSMERROR("type not supported yet");
-        }
-}
-/*}}}*/
-/*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessInternal{{{1*/
-ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessInternal(void){
-
-        /*Nothing added to PVector*/
+        /*No PVector for the Adjoint*/
         return NULL;
-
-}
-/*}}}*/
-/*FUNCTION Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary{{{1*/
-ElementVector* Numericalflux::CreatePVectorAdjointBalancedthicknessBoundary(void){
-
-        /* constants*/
-        const int numdof=NDOF1*NUMVERTICES_BOUNDARY;
-
-        /* Intermediaries*/
-        int        i,j,ig,index1,index2;
-        double     DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN,thickness;
-        double     xyz_list[NUMVERTICES_BOUNDARY][3];
-        double     normal[2];
-        double     L[numdof];
-        GaussTria *gauss;
-
-        /*Initialize Load Vector and return if necessary*/
-        Tria*  tria=(Tria*)element;
-        if(tria->IsOnWater()) return NULL;
-        ElementVector* pe=new ElementVector(nodes,NUMVERTICES_BOUNDARY,this->parameters);
-
-        /*Retrieve all inputs and parameters*/
-        GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES_BOUNDARY);
-        Input* vxaverage_input=tria->inputs->GetInput(VxEnum);           ISSMASSERT(vxaverage_input); 
-        Input* vyaverage_input=tria->inputs->GetInput(VyEnum);           ISSMASSERT(vyaverage_input);
-        Input* thickness_input=tria->inputs->GetInput(ThicknessObsEnum); ISSMASSERT(thickness_input);
-        GetNormal(&normal[0],xyz_list);
-
-        /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/
-        index1=tria->GetNodeIndex(nodes[0]);
-        index2=tria->GetNodeIndex(nodes[1]);
-
-        gauss=new GaussTria();
-        gauss->GaussEdgeCenter(index1,index2);
-        vxaverage_input->GetParameterValue(&mean_vx,gauss);
-        vyaverage_input->GetParameterValue(&mean_vy,gauss);
-        delete gauss;
-        UdotN=mean_vx*normal[0]+mean_vy*normal[1];
-        if (UdotN>0){
-                /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/
-                return NULL;
-        }
-
-        /* Start  looping on the number of gaussian points: */
-        gauss=new GaussTria(index1,index2,2);
-        for(ig=gauss->begin();ig<gauss->end();ig++){
-
-                gauss->GaussPoint(ig);
-
-                tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2);
-
-                vxaverage_input->GetParameterValue(&vx,gauss);
-                vyaverage_input->GetParameterValue(&vy,gauss);
-                thickness_input->GetParameterValue(&thickness,gauss);
-                UdotN=vx*normal[0]+vy*normal[1];
-                tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
-                DL= - gauss->weight*Jdet*UdotN*thickness;
-
-                for(i=0;i<numdof;i++) pe->values[i] += DL*L[i];
-        }
-
-        /*Clean up and return*/
-        delete gauss;
-        return pe;
 }
 /*}}}*/