Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6028) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 6029) @@ -718,6 +718,9 @@ Ke=CreateKMatrixPrognostic(); break; - case BalancedthicknessAnalysisEnum: case AdjointBalancedthicknessAnalysisEnum: + case BalancedthicknessAnalysisEnum: Ke=CreateKMatrixBalancedthickness(); + break; + case AdjointBalancedthicknessAnalysisEnum: + Ke=CreateKMatrixAdjointBalancedthickness(); break; case BalancedvelocitiesAnalysisEnum: @@ -2384,4 +2387,26 @@ /*Tria specific routines: */ +/*FUNCTION Tria::CreateKMatrixAdjointBalancedthickness {{{1*/ +ElementMatrix* Tria::CreateKMatrixAdjointBalancedthickness(void){ + + ElementMatrix* Ke=NULL; + + /*Get Element Matrix of the forward model*/ + switch(GetElementType()){ + case P1Enum: + Ke=CreateKMatrixBalancedthickness_CG(); + break; + case P1DGEnum: + Ke=CreateKMatrixBalancedthickness_DG(); + break; + default: + ISSMERROR("Element type %s not supported yet",EnumToString(GetElementType())); + } + + /*Transpose and return Ke*/ + Ke->Transpose(); + return Ke; +} +/*}}}*/ /*FUNCTION Tria::CreateKMatrixBalancedthickness {{{1*/ ElementMatrix* Tria::CreateKMatrixBalancedthickness(void){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 6028) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 6029) @@ -121,4 +121,5 @@ /*}}}*/ /*Tria specific routines:{{{1*/ + ElementMatrix* CreateKMatrixAdjointBalancedthickness(void); ElementMatrix* CreateKMatrixBalancedthickness(void); ElementMatrix* CreateKMatrixBalancedthickness_DG(void); Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 6028) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 6029) @@ -509,7 +509,7 @@ /*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,6 +532,8 @@ 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,7 +661,7 @@ /*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,6 +683,8 @@ 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,5 +735,7 @@ ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessInternal(void){ - return CreateKMatrixBalancedthicknessInternal(); + ElementMatrix* Ke=CreateKMatrixBalancedthicknessInternal(); + if (Ke) Ke->Transpose(); + return Ke; } /*}}}*/ @@ -737,5 +743,7 @@ ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancedthicknessBoundary(void){ - return CreateKMatrixBalancedthicknessBoundary(); + ElementMatrix* Ke=CreateKMatrixBalancedthicknessBoundary(); + if(Ke) Ke->Transpose(); + return Ke; } /*}}}*/ @@ -779,7 +787,7 @@ /*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,8 +808,11 @@ 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,7 +879,7 @@ /*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,10 +887,5 @@ 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,6 +901,8 @@ 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,87 +933,6 @@ 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();igend();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;ivalues[i] += DL*L[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; } /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 6028) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 6029) @@ -90,6 +90,4 @@ ElementVector* CreatePVectorBalancedthicknessBoundary(void); ElementVector* CreatePVectorAdjointBalancedthickness(void); - ElementVector* CreatePVectorAdjointBalancedthicknessInternal(void); - ElementVector* CreatePVectorAdjointBalancedthicknessBoundary(void); /*}}}*/ Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 6028) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 6029) @@ -43,4 +43,12 @@ this->col_sglobaldoflist=NULL; +} +/*}}}*/ +/*FUNCTION ElementMatrix::ElementMatrix(ElementMatrix* Ke){{{1*/ +ElementMatrix::ElementMatrix(ElementMatrix* Ke){ + + if(!Ke) ISSMERROR("Input Element Matrix is a NULL pointer"); + this->Init(Ke); + return; } /*}}}*/ @@ -319,6 +327,12 @@ /*Intermediaries*/ - double *values_copy=NULL; - int temp; + ElementMatrix* Ke_copy=new ElementMatrix(this); + + /*Update sizes*/ + this->nrows=Ke_copy->ncols; + this->ncols=Ke_copy->nrows; + + /*Transpose values*/ + for (int i=0;inrows;i++) for(int j=0;jncols;j++) this->values[i*this->ncols+j]=Ke_copy->values[j*Ke_copy->ncols+i]; /*Transpose indices*/ @@ -327,16 +341,7 @@ } - /*Transpose values*/ - values_copy=(double*)xmalloc(this->nrows*this->ncols*sizeof(double)); - memcpy(values_copy,this->values,this->nrows*this->ncols*sizeof(double)); - for (int i=0;inrows;i++) for(int j=0;jncols;j++) this->values[j*this->nrows+i]=values_copy[i*this->ncols+j]; - - /*Update sizes*/ - temp=this->nrows; - this->nrows=this->ncols; - this->ncols=temp; - /*Clean up and return*/ - xfree((void**)&values_copy); + delete Ke_copy; + return; } /*}}}*/ Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.h =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.h (revision 6028) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.h (revision 6029) @@ -52,4 +52,5 @@ /*ElementMatrix constructors, destructors {{{1*/ ElementMatrix(); + ElementMatrix(ElementMatrix* Ke); ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2); ElementMatrix(ElementMatrix* Ke1,ElementMatrix* Ke2,ElementMatrix* Ke3);