Index: /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 16861) +++ /issm/trunk-jpl/src/c/analyses/AdjointBalancethicknessAnalysis.cpp (revision 16862) @@ -30,5 +30,102 @@ }/*}}}*/ ElementVector* AdjointBalancethicknessAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + + /*Intermediaries*/ + int meshtype; + Element* basalelement; + + /*Get basal element*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + basalelement = element; + break; + case Mesh3DEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); + break; + default: _error_("mesh "<GetNumberOfNodes(); + + /*Initialize Element vector and vectors*/ + ElementVector* pe = basalelement->NewElementVector(SSAApproximationEnum); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(2*numnodes); + + /*Retrieve all inputs and parameters*/ + basalelement->GetVerticesCoordinates(&xyz_list); + basalelement->FindParam(&num_responses,InversionNumCostFunctionsEnum); + basalelement->FindParam(&responses,NULL,InversionCostFunctionsEnum); + Input* thickness_input = basalelement->GetInput(ThicknessEnum); _assert_(thickness_input); + Input* thicknessobs_input = basalelement->GetInput(InversionThicknessObsEnum); _assert_(thicknessobs_input); + Input* weights_input = basalelement->GetInput(InversionCostFunctionsCoefficientsEnum); _assert_(weights_input); + Input* vx_input = basalelement->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input = basalelement->GetInput(VyEnum); _assert_(vy_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=basalelement->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss); + basalelement->NodalFunctions(basis,gauss); + basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + thickness_input->GetInputValue(&thickness, gauss); + thickness_input->GetInputDerivativeValue(&dH[0],xyz_list,gauss); + thicknessobs_input->GetInputValue(&thicknessobs, gauss); + + /*Loop over all requested responses*/ + for(int resp=0;respGetInputValue(&weight,gauss,responses[resp]); + + switch(responses[resp]){ + case ThicknessAbsMisfitEnum: + for(i=0;ivalues[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i]; + break; + case ThicknessAbsGradientEnum: + for(i=0;ivalues[i]+= - weight*dH[0]*dbasis[0*numnodes+i]*Jdet*gauss->weight; + for(i=0;ivalues[i]+= - weight*dH[1]*dbasis[1*numnodes+i]*Jdet*gauss->weight; + break; + case ThicknessAlongGradientEnum: + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel = sqrt(vx*vx+vy*vy); + vx = vx/(vel+1.e-9); + vy = vy/(vel+1.e-9); + for(i=0;ivalues[i]+= - weight*(dH[0]*vx+dH[1]*vy)*(dbasis[0*numnodes+i]*vx+dbasis[1*numnodes+i]*vy)*Jdet*gauss->weight; + break; + case ThicknessAcrossGradientEnum: + vx_input->GetInputValue(&vx,gauss); + vy_input->GetInputValue(&vy,gauss); + vel = sqrt(vx*vx+vy*vy); + vx = vx/(vel+1.e-9); + vy = vy/(vel+1.e-9); + for(i=0;ivalues[i]+= - weight*(dH[0]*(-vy)+dH[1]*vx)*(dbasis[0*numnodes+i]*(-vy)+dbasis[1*numnodes+i]*vx)*Jdet*gauss->weight; + break; + default: + _error_("response " << EnumToStringx(responses[resp]) << " not supported yet"); + } + } + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + xDelete(dbasis); + if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;}; + delete gauss; + return pe; }/*}}}*/ void AdjointBalancethicknessAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 16861) +++ /issm/trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp (revision 16862) @@ -41,5 +41,5 @@ }/*}}}*/ ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + return NULL; }/*}}}*/ void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16861) +++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp (revision 16862) @@ -121,4 +121,8 @@ case Mesh2DhorizontalEnum: basalelement = element; + break; + case Mesh2DverticalEnum: + if(!element->IsOnBed()) return NULL; + basalelement = element->SpawnBasalElement(); break; case Mesh3DEnum: Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16861) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 16862) @@ -2054,9 +2054,12 @@ /*Retrieve all inputs and parameters*/ element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); Input* loadingforcex_input=element->GetInput(LoadingforceXEnum); _assert_(loadingforcex_input); Input* loadingforcey_input=element->GetInput(LoadingforceYEnum); _assert_(loadingforcey_input); - Input* loadingforcez_input=element->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); - IssmDouble rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble gravity =element->GetMaterialParameter(ConstantsGEnum); + Input* loadingforcez_input=NULL; + if(dim==3){ + loadingforcez_input=element->GetInput(LoadingforceZEnum); _assert_(loadingforcez_input); + } /* Start looping on the number of gaussian points: */ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16861) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 16862) @@ -514,2 +514,66 @@ } /*}}}*/ +/*FUNCTION Seg::DeleteMaterials{{{*/ +void Seg::DeleteMaterials(void){ + delete this->material; +} +/*}}}*/ +/*FUNCTION Seg::GetInput(int inputenum) {{{*/ +Input* Seg::GetInput(int inputenum){ + return inputs->GetInput(inputenum); +} +/*}}}*/ +/*FUNCTION Seg::GetNumberOfNodes;{{{*/ +int Seg::GetNumberOfNodes(void){ + return this->NumberofNodes(); +} +/*}}}*/ +/*FUNCTION Seg::GetVerticesCoordinates(IssmDouble** pxyz_list){{{*/ +void Seg::GetVerticesCoordinates(IssmDouble** pxyz_list){ + + IssmDouble* xyz_list = xNew(NUMVERTICES*3); + ::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES); + + /*Assign output pointer*/ + *pxyz_list = xyz_list; + +}/*}}}*/ +/*FUNCTION Seg::JacobianDeterminant{{{*/ +void Seg::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussSegEnum); + this->GetJacobianDeterminant(pJdet,xyz_list,(GaussSeg*)gauss); + +} +/*}}}*/ +/*FUNCTION Seg::NewGauss(int order){{{*/ +Gauss* Seg::NewGauss(int order){ + return new GaussSeg(order); +} +/*}}}*/ +/*FUNCTION Seg::NewElementVector{{{*/ +ElementVector* Seg::NewElementVector(int approximation_enum){ + return new ElementVector(nodes,this->NumberofNodes(),this->parameters,approximation_enum); +} +/*}}}*/ +/*FUNCTION Seg::NewElementMatrix{{{*/ +ElementMatrix* Seg::NewElementMatrix(int approximation_enum){ + return new ElementMatrix(nodes,this->NumberofNodes(),this->parameters,approximation_enum); +} +/*}}}*/ +/*FUNCTION Seg::NodalFunctions{{{*/ +void Seg::NodalFunctions(IssmDouble* basis, Gauss* gauss){ + + _assert_(gauss->Enum()==GaussSegEnum); + this->GetNodalFunctions(basis,(GaussSeg*)gauss); + +} +/*}}}*/ +/*FUNCTION Seg::NodalFunctionsDerivatives{{{*/ +void Seg::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){ + + _assert_(gauss->Enum()==GaussSegEnum); + this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussSeg*)gauss); + +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16861) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16862) @@ -77,5 +77,5 @@ void ComputeStressTensor(){_error_("not implemented yet");}; void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; - void DeleteMaterials(void){_error_("not implemented yet");}; + void DeleteMaterials(void); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");}; @@ -104,9 +104,9 @@ void GetNodesSidList(int* sidlist){_error_("not implemented yet");}; void GetNodesLidList(int* lidlist){_error_("not implemented yet");}; - int GetNumberOfNodes(void){_error_("not implemented yet");}; + int GetNumberOfNodes(void); int GetNumberOfNodesVelocity(void){_error_("not implemented yet");}; int GetNumberOfNodesPressure(void){_error_("not implemented yet");}; int GetNumberOfVertices(void){_error_("not implemented yet");}; - void GetVerticesCoordinates(IssmDouble** pxyz_list){_error_("not implemented yet");}; + void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");}; void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; @@ -117,5 +117,5 @@ bool IsFloating(){_error_("not implemented yet");}; bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; - void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void JacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss); void JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; @@ -123,8 +123,8 @@ void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){_error_("not implemented yet");}; IssmDouble MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");}; - void NodalFunctions(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctions(IssmDouble* basis,Gauss* gauss); void NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");}; - void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; + void NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; bool NoIceInElement(){_error_("not implemented yet");}; @@ -141,5 +141,5 @@ void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; int VelocityInterpolation(void){_error_("not implemented yet");}; - Input* GetInput(int inputenum){_error_("not implemented yet");}; + Input* GetInput(int inputenum); Input* GetMaterialInput(int inputenum){_error_("not implemented yet");}; void GetInputListOnVertices(IssmDouble* pvalue,int enumtype){_error_("not implemented yet");}; @@ -152,5 +152,5 @@ IssmDouble GetZcoord(Gauss* gauss){_error_("not implemented yet");}; Gauss* NewGauss(void){_error_("not implemented yet");}; - Gauss* NewGauss(int order){_error_("not implemented yet");}; + Gauss* NewGauss(int order); Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");}; @@ -158,6 +158,6 @@ Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");}; Gauss* NewGaussTop(int order){_error_("not implemented yet");}; - ElementVector* NewElementVector(int approximation_enum){_error_("not implemented yet");}; - ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");}; + ElementVector* NewElementVector(int approximation_enum); + ElementMatrix* NewElementMatrix(int approximation_enum); int VertexConnectivity(int vertexindex){_error_("not implemented yet");}; void VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16861) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 16862) @@ -2249,4 +2249,12 @@ } /*}}}*/ +/*FUNCTION Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){{{*/ +Gauss* Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){ + + IssmDouble area_coordinates[4][3]; + GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4); + return new GaussTria(area_coordinates,order_vert); +} +/*}}}*/ /*FUNCTION Tria::NewElementVector{{{*/ ElementVector* Tria::NewElementVector(int approximation_enum){ @@ -2677,13 +2685,18 @@ Element* Tria::SpawnBasalElement(void){ + int index1,index2; int meshtype; - /*go into parameters and get the analysis_counter: */ this->parameters->FindParam(&meshtype,MeshTypeEnum); - - if(meshtype==Mesh2DhorizontalEnum){ - return this; - } - else _error_("not implemented yet"); + switch(meshtype){ + case Mesh2DhorizontalEnum: + return this; + case Mesh2DverticalEnum: + _assert_(HasEdgeOnBed()); + this->EdgeOnBedIndices(&index1,&index2); + return SpawnSeg(index1,index2); + default: + _error_("not implemented yet"); + } } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16861) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16862) @@ -296,5 +296,5 @@ Gauss* NewGauss(int order); Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order); - Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");}; + Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert); Gauss* NewGaussBase(int order){_error_("not implemented yet");}; Gauss* NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};