Index: /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17516) @@ -214,6 +214,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numdof = vnumnodes*3 + pnumnodes; @@ -320,6 +320,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ @@ -939,6 +939,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int vnumdof = vnumnodes*3; int pnumdof = pnumnodes*1; Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17516) @@ -73,5 +73,5 @@ void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ - iomodel->FetchData(1,MeshVertexonbedEnum); + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbedEnum); ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum); iomodel->DeleteData(1,MeshVertexonbedEnum); @@ -87,18 +87,5 @@ void DamageEvolutionAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - /*create penalties for nodes: no node can have a damage > 1*/ - iomodel->FetchData(1,DamageSpcdamageEnum); - CreateSingleNodeToElementConnectivity(iomodel); - - for(int i=0;inumberofvertices;i++){ - - /*keep only this partition's nodes:*/ - if(iomodel->my_vertices[i]){ - if (xIsNan(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes! - loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum)); - } - } - } - iomodel->DeleteData(1,DamageSpcdamageEnum); + /*Nothing for now*/ }/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 17516) @@ -60,5 +60,5 @@ IssmDouble *nodeonbed=NULL; iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); - iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); for(int i=0;i(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + if(iomodel->meshtype!=Mesh2DhorizontalEnum){ + if(!(reCast(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + } /*Get node ids*/ Index: /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 17516) @@ -35,5 +35,5 @@ iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); iomodel->FetchDataToInput(elements,VxEnum); - iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); if(iomodel->meshtype==Mesh3DEnum){ iomodel->FetchDataToInput(elements,VzEnum); @@ -69,5 +69,5 @@ IssmDouble *nodeonsurface=NULL; iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); - iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); for(int i=0;i(nodeonsurface[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonsurface[reCast(vertex_pairing[2*i+1])-1]))) continue; + if(iomodel->meshtype!=Mesh2DhorizontalEnum){ + if(!(reCast(nodeonsurface[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonsurface[reCast(vertex_pairing[2*i+1])-1]))) continue; + } /*Get node ids*/ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17516) @@ -195,5 +195,5 @@ IssmDouble *nodeonbed=NULL; iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); - iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); + if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); for(int i=0;i(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + if(iomodel->meshtype!=Mesh2DhorizontalEnum){ + if(!(reCast(nodeonbed[reCast(vertex_pairing[2*i+0])-1])) || !(reCast(nodeonbed[reCast(vertex_pairing[2*i+1])-1]))) continue; + } /*Get node ids*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17516) @@ -2754,6 +2754,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Initialize output vector*/ @@ -2784,6 +2784,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; @@ -2886,6 +2886,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numdof = vnumnodes*dim + pnumnodes; int bsize = epssize + 2; @@ -2970,6 +2970,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numdof = vnumnodes*dim + pnumnodes; @@ -3065,6 +3065,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numdof = vnumnodes*dim + pnumnodes; @@ -3126,6 +3126,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ @@ -3209,6 +3209,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ @@ -3288,6 +3288,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ @@ -3365,6 +3365,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numvertices = element->GetNumberOfVertices(); @@ -3659,6 +3659,6 @@ /*Fetch number of nodes for this finite element*/ - int pnumnodes = element->GetNumberOfNodesPressure(); - int vnumnodes = element->GetNumberOfNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); int pnumdof = pnumnodes; int vnumdof = vnumnodes*dim; @@ -3786,6 +3786,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int vnumdof = vnumnodes*dim; int pnumdof = pnumnodes*1; @@ -4321,6 +4321,6 @@ if(element->IsFloating() || !element->IsOnBed()) return NULL; - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numnodes = 2*vnumnodes-1+pnumnodes; @@ -4439,6 +4439,6 @@ Element* basaltria=pentabase->SpawnBasalElement(); - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numnodes = 2*vnumnodes-1+pnumnodes; @@ -4608,6 +4608,6 @@ if(approximation!=HOFSApproximationEnum) return NULL; - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); int numnodes = vnumnodes+pnumnodes; @@ -4687,6 +4687,6 @@ element->GetInputValue(&approximation,ApproximationEnum); if(approximation!=HOFSApproximationEnum) return NULL; - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ @@ -4771,6 +4771,6 @@ element->GetInputValue(&approximation,ApproximationEnum); if(approximation!=SSAFSApproximationEnum) return NULL; - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ @@ -4847,6 +4847,6 @@ element->GetInputValue(&approximation,ApproximationEnum); if(approximation!=SSAFSApproximationEnum) return NULL; - int vnumnodes = element->GetNumberOfNodesVelocity(); - int pnumnodes = element->GetNumberOfNodesPressure(); + int vnumnodes = element->NumberofNodesVelocity(); + int pnumnodes = element->NumberofNodesPressure(); /*Prepare coordinate system list*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17516) @@ -17,4 +17,6 @@ /*Constructors/destructor/copy*/ Element::Element(){/*{{{*/ + this->id = -1; + this->sid = -1; this->inputs = NULL; this->nodes = NULL; @@ -130,4 +132,69 @@ delete this->material; }/*}}}*/ +void Element::DeepEcho(void){/*{{{*/ + + _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n"); + _printf_(" id : "<id <<"\n"); + _printf_(" sid: "<sid<<"\n"); + if(vertices){ + int numvertices = this->GetNumberOfVertices(); + for(int i=0;iEcho(); + } + else _printf_("vertices = NULL\n"); + + if(nodes){ + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iDeepEcho(); + } + else _printf_("nodes = NULL\n"); + + if (material) material->DeepEcho(); + else _printf_("material = NULL\n"); + + if (matpar) matpar->DeepEcho(); + else _printf_("matpar = NULL\n"); + + _printf_(" parameters\n"); + if (parameters) parameters->DeepEcho(); + else _printf_("parameters = NULL\n"); + + _printf_(" inputs\n"); + if (inputs) inputs->DeepEcho(); + else _printf_("inputs=NULL\n"); + + return; +} +/*}}}*/ +void Element::Echo(void){/*{{{*/ + _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n"); + _printf_(" id : "<id <<"\n"); + _printf_(" sid: "<sid<<"\n"); + if(vertices){ + int numvertices = this->GetNumberOfVertices(); + for(int i=0;iEcho(); + } + else _printf_("vertices = NULL\n"); + + if(nodes){ + int numnodes = this->GetNumberOfNodes(); + for(int i=0;iEcho(); + } + else _printf_("nodes = NULL\n"); + + if (material) material->Echo(); + else _printf_("material = NULL\n"); + + if (matpar) matpar->Echo(); + else _printf_("matpar = NULL\n"); + + _printf_(" parameters\n"); + if (parameters) parameters->Echo(); + else _printf_("parameters = NULL\n"); + + _printf_(" inputs\n"); + if (inputs) inputs->Echo(); + else _printf_("inputs=NULL\n"); +} +/*}}}*/ IssmDouble Element::Divergence(void){/*{{{*/ /*Compute element divergence*/ @@ -162,4 +229,16 @@ return divergence; }/*}}}*/ +void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ + matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); +}/*}}}*/ +void Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ + matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure); +}/*}}}*/ +IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ + return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); +}/*}}}*/ +IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/ + return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); +}/*}}}*/ void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/ this->parameters->FindParam(pvalue,paramenum); @@ -200,5 +279,5 @@ /*Fetch number of nodes and dof for this finite element*/ - int numnodes = this->GetNumberOfNodesVelocity(); + int numnodes = this->NumberofNodesVelocity(); /*First, figure out size of doflist and create it: */ @@ -223,6 +302,6 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->GetNumberOfNodesVelocity(); - int pnumnodes = this->GetNumberOfNodesPressure(); + int vnumnodes = this->NumberofNodesVelocity(); + int pnumnodes = this->NumberofNodesPressure(); /*First, figure out size of doflist and create it: */ @@ -388,5 +467,5 @@ _assert_(pvalue); - int numnodes = this->GetNumberOfNodesVelocity(); + int numnodes = this->NumberofNodesVelocity(); Input *input = this->GetInput(enumtype); if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); @@ -468,4 +547,16 @@ int numvertices = this->GetNumberOfVertices(); for(int i=0;ivertices[i]->Connectivity(); +} +/*}}}*/ +bool Element::HasNodeOnBed(){/*{{{*/ + return (this->inputs->Max(MeshVertexonbedEnum)>0.); +}/*}}}*/ +bool Element::HasNodeOnSurface(){/*{{{*/ + return (this->inputs->Max(MeshVertexonsurfaceEnum)>0.); +}/*}}}*/ +int Element::Id(){/*{{{*/ + + return this->id; + } /*}}}*/ @@ -587,4 +678,8 @@ return shelf; }/*}}}*/ +bool Element::IsIceInElement(){/*{{{*/ + return (this->inputs->Min(MaskIceLevelsetEnum)<0.); +} +/*}}}*/ bool Element::IsInput(int name){/*{{{*/ if ( @@ -675,4 +770,7 @@ } /*}}}*/ +IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/ + return this->matpar->PureIceEnthalpy(pressure); +}/*}}}*/ void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/ @@ -773,4 +871,80 @@ } /*}}}*/ +void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/ + + /*Intermediaries*/ + const int numnodes = this->GetNumberOfNodes(); + + /*Output */ + int d_nz = 0; + int o_nz = 0; + + /*Loop over all nodes*/ + for(int i=0;inodes[i]->Lid()]){ + + /*flag current node so that no other element processes it*/ + flags[this->nodes[i]->Lid()]=true; + + int counter=0; + while(flagsindices[counter]>=0) counter++; + flagsindices[counter]=this->nodes[i]->Lid(); + + /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ + switch(set2_enum){ + case FsetEnum: + if(nodes[i]->indexing.fsize){ + if(this->nodes[i]->IsClone()) + o_nz += 1; + else + d_nz += 1; + } + break; + case GsetEnum: + if(nodes[i]->indexing.gsize){ + if(this->nodes[i]->IsClone()) + o_nz += 1; + else + d_nz += 1; + } + break; + case SsetEnum: + if(nodes[i]->indexing.ssize){ + if(this->nodes[i]->IsClone()) + o_nz += 1; + else + d_nz += 1; + } + break; + default: _error_("not supported"); + } + } + } + + /*Special case: 2d/3d coupling, the node of this element might be connected + *to the basal element*/ + int analysis_type,approximation,numlayers; + parameters->FindParam(&analysis_type,AnalysisTypeEnum); + if(analysis_type==StressbalanceAnalysisEnum){ + inputs->GetInputValue(&approximation,ApproximationEnum); + if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){ + parameters->FindParam(&numlayers,MeshNumberoflayersEnum); + o_nz += numlayers*3; + d_nz += numlayers*3; + } + } + + /*Assign output pointers: */ + *pd_nz=d_nz; + *po_nz=o_nz; +} +/*}}}*/ +int Element::Sid(){/*{{{*/ + + return this->sid; + +} +/*}}}*/ IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/ _assert_(matpar); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17516) @@ -38,4 +38,6 @@ public: + int id; + int sid; Inputs *inputs; Node **nodes; @@ -55,6 +57,12 @@ /* bool AnyActive(void); */ void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array); + void Echo(); + void DeepEcho(); void DeleteMaterials(void); IssmDouble Divergence(void); + void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); + void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); + IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); + IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); void FindParam(bool* pvalue,int paramenum); void FindParam(int* pvalue,int paramenum); @@ -81,4 +89,8 @@ void GetVerticesSidList(int* sidlist); void GetVerticesConnectivityList(int* connectivitylist); + bool HasNodeOnBed(); + bool HasNodeOnSurface(); + int Id(); + int Sid(); void InputChangeName(int enum_type,int enum_type_old); void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); @@ -86,4 +98,5 @@ void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); + bool IsIceInElement(); bool IsInput(int name); bool IsFloating(); @@ -91,7 +104,9 @@ ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); + IssmDouble PureIceEnthalpy(IssmDouble pressure); void ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum); void ResultToVector(Vector* vector,int output_enum); void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum); + void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); void StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input); @@ -135,10 +150,6 @@ virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0; - virtual void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; - virtual void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure)=0; - virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0; - virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0; - virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0; + virtual int FiniteElement(void)=0; virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0; @@ -154,5 +165,4 @@ virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0; virtual void NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0; - virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0; virtual Element* GetUpperElement(void)=0; @@ -169,14 +179,9 @@ virtual int GetNodeIndex(Node* node)=0; virtual int GetNumberOfNodes(void)=0; - virtual int GetNumberOfNodesVelocity(void)=0; - virtual int GetNumberOfNodesPressure(void)=0; virtual int GetNumberOfVertices(void)=0; - virtual int Id()=0; - virtual int Sid()=0; virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; virtual bool IsOnBed()=0; virtual bool IsOnSurface()=0; - virtual bool IsIceInElement()=0; virtual void GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0; virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17516) @@ -371,24 +371,4 @@ } /*}}}*/ -/*FUNCTION Penta::DeepEcho{{{*/ -void Penta::DeepEcho(void){ - - _printf_("Penta:\n"); - _printf_(" id: " << id << "\n"); - nodes[0]->DeepEcho(); - nodes[1]->DeepEcho(); - nodes[2]->DeepEcho(); - nodes[3]->DeepEcho(); - nodes[4]->DeepEcho(); - nodes[5]->DeepEcho(); - material->DeepEcho(); - matpar->DeepEcho(); - _printf_(" neighbor ids: " << verticalneighbors[0]->Id() << "-" << verticalneighbors[1]->Id() << "\n"); - _printf_(" parameters\n"); - parameters->DeepEcho(); - _printf_(" inputs\n"); - inputs->DeepEcho(); -} -/*}}}*/ /*FUNCTION Penta::Delta18oParameterization{{{*/ void Penta::Delta18oParameterization(void){ @@ -463,30 +443,4 @@ /*clean-up*/ delete gauss; -} -/*}}}*/ -/*FUNCTION Penta::Echo{{{*/ - -void Penta::Echo(void){ - this->DeepEcho(); -} -/*}}}*/ -/*FUNCTION Penta::ThermalToEnthalpy{{{*/ -void Penta::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){ - matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); -} -/*}}}*/ -/*FUNCTION Penta::EnthalpyToThermal{{{*/ -void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){ - matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure); -} -/*}}}*/ -/*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/ -IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){ - return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); -} -/*}}}*/ -/*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/ -IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){ - return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); } /*}}}*/ @@ -848,14 +802,4 @@ } /*}}}*/ -/*FUNCTION Penta::GetNumberOfNodesPressure;{{{*/ -int Penta::GetNumberOfNodesPressure(void){ - return this->NumberofNodesPressure(); -} -/*}}}*/ -/*FUNCTION Penta::GetNumberOfNodesVelocity;{{{*/ -int Penta::GetNumberOfNodesVelocity(void){ - return this->NumberofNodesVelocity(); -} -/*}}}*/ /*FUNCTION Penta::GetNumberOfVertices;{{{*/ int Penta::GetNumberOfVertices(void){ @@ -1160,16 +1104,4 @@ /*Assign output pointer*/ *pxyz_zero= xyz_zero; -} -/*}}}*/ -/*FUNCTION Penta::Sid {{{*/ -int Penta::Sid(){ - - return sid; - -} -/*}}}*/ -/*FUNCTION Penta::Id {{{*/ -int Penta::Id(void){ - return id; } /*}}}*/ @@ -1643,18 +1575,4 @@ } /*}}}*/ -/*FUNCTION Penta::IsIceInElement {{{*/ -bool Penta::IsIceInElement(){ - - /*Get levelset*/ - IssmDouble ls[NUMVERTICES]; - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /*If the level set one one of the nodes is <0, ice is present in this element*/ - if(ls[0]<0. || ls[1]<0. || ls[2]<0.) - return true; - else - return false; -} -/*}}}*/ /*FUNCTION Penta::MinEdgeLength{{{*/ IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){ @@ -1813,10 +1731,9 @@ void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ - int i; IssmDouble v13[3],v23[3]; IssmDouble normal[3]; IssmDouble normal_norm; - for (i=0;i<3;i++){ + for(int i=0;i<3;i++){ v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; @@ -1826,5 +1743,5 @@ normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; - normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2) ); + normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]); /*Bed normal is opposite to surface normal*/ @@ -1917,10 +1834,4 @@ /*clean-up*/ delete gauss; -} -/*}}}*/ -/*FUNCTION Penta::PureIceEnthalpy{{{*/ -IssmDouble Penta::PureIceEnthalpy(IssmDouble pressure){ - - return this->matpar->PureIceEnthalpy(pressure); } /*}}}*/ @@ -2070,75 +1981,4 @@ void Penta::SetTemporaryElementType(int element_type_in){ this->element_type=element_type_in; -} -/*}}}*/ -/*FUNCTION Penta::SetwiseNodeConnectivity{{{*/ -void Penta::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){ - - /*Intermediaries*/ - const int numnodes = this->NumberofNodes(); - - /*Output */ - int d_nz = 0; - int o_nz = 0; - - /*Loop over all nodes*/ - for(int i=0;inodes[i]->Lid()]){ - - /*flag current node so that no other element processes it*/ - flags[this->nodes[i]->Lid()]=true; - - int counter=0; - while(flagsindices[counter]>=0) counter++; - flagsindices[counter]=this->nodes[i]->Lid(); - - /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ - switch(set2_enum){ - case FsetEnum: - if(nodes[i]->indexing.fsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case GsetEnum: - if(nodes[i]->indexing.gsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case SsetEnum: - if(nodes[i]->indexing.ssize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - default: _error_("not supported"); - } - } - } - - /*Special case: 2d/3d coupling, the node of this element might be connected - *to the basal element*/ - int analysis_type,approximation,numlayers; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - if(analysis_type==StressbalanceAnalysisEnum){ - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){ - parameters->FindParam(&numlayers,MeshNumberoflayersEnum); - o_nz += numlayers*3; - d_nz += numlayers*3; - } - } - - /*Assign output pointers: */ - *pd_nz=d_nz; - *po_nz=o_nz; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17516) @@ -32,7 +32,4 @@ public: - int id; - int sid; - Penta **verticalneighbors; // 2 neighbors: first one under, second one above @@ -44,8 +41,5 @@ /*Object virtual functions definitions: {{{*/ Object *copy(); - void DeepEcho(); - void Echo(); int ObjectEnum(); - int Id(); /*}}}*/ /*Update virtual functions definitions: {{{*/ @@ -67,8 +61,5 @@ int FiniteElement(void); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); - void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void Delta18oParameterization(void); - void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); - void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); Penta* GetUpperPenta(void); Penta* GetLowerPenta(void); @@ -83,6 +74,4 @@ int GetNodeIndex(Node* node); int GetNumberOfNodes(void); - int GetNumberOfNodesPressure(void); - int GetNumberOfNodesVelocity(void); int GetNumberOfVertices(void); void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); @@ -94,5 +83,4 @@ void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); - int Sid(); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputDuplicate(int original_enum,int new_enum); @@ -101,5 +89,4 @@ int NumberofNodesPressure(void); int VelocityInterpolation(); - IssmDouble PureIceEnthalpy(IssmDouble pressure); int PressureInterpolation(); bool IsZeroLevelset(int levelset_enum); @@ -189,6 +176,4 @@ void NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list); ElementMatrix* CreateBasalMassMatrix(void); - IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); - IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); @@ -208,5 +193,4 @@ void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); - bool IsIceInElement(void); Gauss* NewGauss(void); Gauss* NewGauss(int order); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17516) @@ -65,29 +65,4 @@ } /*}}}*/ -/*FUNCTION Seg::Echo{{{*/ -void Seg::Echo(void){ - _printf_("Seg:\n"); - _printf_(" id: " << id << "\n"); - if(nodes){ - nodes[0]->Echo(); - nodes[1]->Echo(); - } - else _printf_("nodes = NULL\n"); - - if (material) material->Echo(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->Echo(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->Echo(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->Echo(); - else _printf_("inputs=NULL\n"); -} -/*}}}*/ /*FUNCTION Seg::FiniteElement{{{*/ int Seg::FiniteElement(void){ @@ -95,43 +70,8 @@ } /*}}}*/ -/*FUNCTION Seg::DeepEcho{{{*/ -void Seg::DeepEcho(void){ - - _printf_("Seg:\n"); - _printf_(" id: " << id << "\n"); - if(nodes){ - nodes[0]->DeepEcho(); - nodes[1]->DeepEcho(); - } - else _printf_("nodes = NULL\n"); - - if (material) material->DeepEcho(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->DeepEcho(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->DeepEcho(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->DeepEcho(); - else _printf_("inputs=NULL\n"); - - return; -} -/*}}}*/ /*FUNCTION Seg::ObjectEnum{{{*/ int Seg::ObjectEnum(void){ return SegEnum; - -} -/*}}}*/ -/*FUNCTION Seg::Id {{{*/ -int Seg::Id(){ - - return id; } @@ -186,18 +126,4 @@ }/*}}}*/ -/*FUNCTION Seg::IsIceInElement {{{*/ -bool Seg::IsIceInElement(){ - - /*Get levelset*/ - IssmDouble ls[NUMVERTICES]; - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /*If the level set on one of the nodes is <0, ice is present in this element*/ - if(ls[0]<0. || ls[1]<0.) - return true; - else - return false; -} -/*}}}*/ bool Seg::IsIcefront(void){/*{{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17516) @@ -30,7 +30,4 @@ public: - int id; - int sid; - /*Seg constructors, destructors {{{*/ Seg(){}; @@ -39,7 +36,4 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); int ObjectEnum(); Object *copy(); @@ -64,11 +58,6 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; 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");}; void Delta18oParameterization(void){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; - void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");}; - void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; - IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; - IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; int FiniteElement(void); Element* GetUpperElement(void){_error_("not implemented yet");}; @@ -78,11 +67,8 @@ int GetNodeIndex(Node* node){_error_("not implemented yet");}; int GetNumberOfNodes(void); - int GetNumberOfNodesVelocity(void){_error_("not implemented yet");}; - int GetNumberOfNodesPressure(void){_error_("not implemented yet");}; int GetNumberOfVertices(void); void GetVerticesCoordinates(IssmDouble** pxyz_list); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");}; void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; - int Sid(){_error_("not implemented yet");}; bool IsOnBed(){_error_("not implemented yet");}; bool IsOnSurface(){_error_("not implemented yet");}; @@ -102,5 +88,4 @@ void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; - bool IsIceInElement(); void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; @@ -111,5 +96,4 @@ Element* SpawnTopElement(void){_error_("not implemented yet");}; IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; - IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");}; int PressureInterpolation(void){_error_("not implemented yet");}; void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17516) @@ -52,66 +52,8 @@ /*}}}*/ -/*FUNCTION Tetra::Echo{{{*/ -void Tetra::Echo(void){ - _printf_("Tetra:\n"); - _printf_(" id: " << id << "\n"); - if(nodes){ - nodes[0]->Echo(); - nodes[1]->Echo(); - nodes[2]->Echo(); - nodes[3]->Echo(); - } - else _printf_("nodes = NULL\n"); - - if (material) material->Echo(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->Echo(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->Echo(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->Echo(); - else _printf_("inputs=NULL\n"); -} -/*}}}*/ /*FUNCTION Tetra::FiniteElement{{{*/ int Tetra::FiniteElement(void){ return this->element_type; -} -/*}}}*/ -/*FUNCTION Tetra::DeepEcho{{{*/ -void Tetra::DeepEcho(void){ - - _printf_("Tetra:\n"); - _printf_(" id: " << id << "\n"); - if(nodes){ - nodes[0]->DeepEcho(); - nodes[1]->DeepEcho(); - nodes[2]->DeepEcho(); - nodes[3]->DeepEcho(); - } - else _printf_("nodes = NULL\n"); - - if (material) material->DeepEcho(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->DeepEcho(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->DeepEcho(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->DeepEcho(); - else _printf_("inputs=NULL\n"); - - return; -} -/*}}}*/ +} /*}}}*/ /*FUNCTION Tetra::ObjectEnum{{{*/ int Tetra::ObjectEnum(void){ @@ -119,11 +61,5 @@ return TetraEnum; -} -/*}}}*/ -/*FUNCTION Tetra::Id {{{*/ -int Tetra::Id(){ - return id; -} -/*}}}*/ +}/*}}}*/ /*FUNCTION Tetra::AddInput{{{*/ @@ -240,14 +176,4 @@ } /*}}}*/ -/*FUNCTION Tetra::GetNumberOfNodesPressure THIS ONE (and corresponding TetraRef function){{{*/ -int Tetra::GetNumberOfNodesPressure(void){ - return this->NumberofNodesPressure(); -} -/*}}}*/ -/*FUNCTION Tetra::GetNumberOfNodesVelocity;{{{*/ -int Tetra::GetNumberOfNodesVelocity(void){ - return this->NumberofNodesVelocity(); -} -/*}}}*/ /*FUNCTION Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/ void Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){ @@ -339,17 +265,4 @@ return false; } -} -/*}}}*/ -/*FUNCTION Tetra::HasNodeOnBed THIS ONE{{{*/ -bool Tetra::HasNodeOnBed(){ - - IssmDouble values[NUMVERTICES]; - IssmDouble sum; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&values[0],MeshVertexonbedEnum); - sum = values[0]+values[1]+values[2]+values[3]; - - return sum>0.; } /*}}}*/ @@ -485,19 +398,4 @@ } /*}}}*/ -/*FUNCTION Tetra::IsIceInElement THIS ONE{{{*/ -bool Tetra::IsIceInElement(){ - - /*Get levelset*/ - IssmDouble ls[NUMVERTICES]; - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /*If the level set on one of the nodes is <0, ice is present in this element*/ - for(int i=0;ielement_type==MINIcondensedEnum){ - _error_("Not implemented"); + int indices[3]={12,13,14}; + pe->StaticCondensation(Ke,3,&indices[0]); } else if(this->element_type==P1bubblecondensedEnum){ @@ -681,5 +625,6 @@ if(Ke){ if(this->element_type==MINIcondensedEnum){ - _error_("Not implemented"); + int indices[3]={12,13,14}; + Ke->StaticCondensation(3,&indices[0]); } else if(this->element_type==P1bubblecondensedEnum){ @@ -768,69 +713,4 @@ if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); else this->nodes=NULL; - -} -/*}}}*/ -/*FUNCTION Tetra::SetwiseNodeConnectivity THIS ONE (take from Penta.cpp){{{*/ -void Tetra::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){ - - /*Intermediaries*/ - const int numnodes = this->NumberofNodes(); - - /*Output */ - int d_nz = 0; - int o_nz = 0; - - /*Loop over all nodes*/ - for(int i=0;inodes[i]->Lid()]){ - - /*flag current node so that no other element processes it*/ - flags[this->nodes[i]->Lid()]=true; - - int counter=0; - while(flagsindices[counter]>=0) counter++; - flagsindices[counter]=this->nodes[i]->Lid(); - - /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ - switch(set2_enum){ - case FsetEnum: - if(nodes[i]->indexing.fsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case GsetEnum: - if(nodes[i]->indexing.gsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case SsetEnum: - if(nodes[i]->indexing.ssize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - default: _error_("not supported"); - } - } - } - - /*Assign output pointers: */ - *pd_nz=d_nz; - *po_nz=o_nz; -} -/*}}}*/ -/*FUNCTION Tetra::Sid THIS ONE{{{*/ -int Tetra::Sid(){ - - return sid; } Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17516) @@ -30,7 +30,4 @@ public: - int id; - int sid; - /*Tetra constructors, destructors {{{*/ Tetra(){}; @@ -39,7 +36,4 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); int ObjectEnum(); Object *copy(); @@ -64,11 +58,6 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); - void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void Delta18oParameterization(void){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; - void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");}; - void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; - IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; - IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; void FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3); void FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3); @@ -81,6 +70,4 @@ int GetNodeIndex(Node* node){_error_("not implemented yet");}; int GetNumberOfNodes(void); - int GetNumberOfNodesVelocity(void); - int GetNumberOfNodesPressure(void); int GetNumberOfVertices(void); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); @@ -88,6 +75,4 @@ bool HasFaceOnBed(); bool HasFaceOnSurface(); - bool HasNodeOnBed(); - int Sid(); bool IsOnBed(); bool IsOnSurface(); @@ -107,8 +92,7 @@ void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); - bool IsIceInElement(); void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); - void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; - void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; + void NormalTop(IssmDouble* normal,IssmDouble* xyz_list); + void NormalBase(IssmDouble* normal,IssmDouble* xyz_list); int NumberofNodesVelocity(void); int NumberofNodesPressure(void); @@ -117,5 +101,4 @@ Tria* SpawnTria(int index1,int index2,int index3); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; - IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");}; int PressureInterpolation(void); void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");}; @@ -228,3 +211,3 @@ /*}}}*/ }; -#endif /* _SEG_H */ +#endif /* _TETRA_H_*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17516) @@ -95,62 +95,4 @@ /*Other*/ -/*FUNCTION Tria::SetwiseNodeConnectivity{{{*/ -void Tria::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){ - - /*Intermediaries*/ - const int numnodes = this->NumberofNodes(); - - /*Output */ - int d_nz = 0; - int o_nz = 0; - - /*Loop over all nodes*/ - for(int i=0;inodes[i]->Lid()]){ - - /*flag current node so that no other element processes it*/ - flags[this->nodes[i]->Lid()]=true; - - int counter=0; - while(flagsindices[counter]>=0) counter++; - flagsindices[counter]=this->nodes[i]->Lid(); - - /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ - switch(set2_enum){ - case FsetEnum: - if(nodes[i]->indexing.fsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case GsetEnum: - if(nodes[i]->indexing.gsize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - case SsetEnum: - if(nodes[i]->indexing.ssize){ - if(this->nodes[i]->IsClone()) - o_nz += 1; - else - d_nz += 1; - } - break; - default: _error_("not supported"); - } - } - } - - /*Assign output pointers: */ - *pd_nz=d_nz; - *po_nz=o_nz; -} -/*}}}*/ /*FUNCTION Tria::AddBasalInput{{{*/ void Tria::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){ @@ -292,33 +234,4 @@ } /*}}}*/ -/*FUNCTION Tria::DeepEcho{{{*/ -void Tria::DeepEcho(void){ - - _printf_("Tria:\n"); - _printf_(" id: " << id << "\n"); - if(nodes){ - nodes[0]->DeepEcho(); - nodes[1]->DeepEcho(); - nodes[2]->DeepEcho(); - } - else _printf_("nodes = NULL\n"); - - if (material) material->DeepEcho(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->DeepEcho(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->DeepEcho(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->DeepEcho(); - else _printf_("inputs=NULL\n"); - - return; -} -/*}}}*/ /*FUNCTION Tria::Delta18oParameterization{{{*/ void Tria::Delta18oParameterization(void){ @@ -414,30 +327,4 @@ *hy=ymax-ymin; *hz=0.; -} -/*}}}*/ -/*FUNCTION Tria::Echo{{{*/ -void Tria::Echo(void){ - _printf_("Tria:\n"); - _printf_(" id: " << id << "\n"); - if(nodes){ - nodes[0]->Echo(); - nodes[1]->Echo(); - nodes[2]->Echo(); - } - else _printf_("nodes = NULL\n"); - - if (material) material->Echo(); - else _printf_("material = NULL\n"); - - if (matpar) matpar->Echo(); - else _printf_("matpar = NULL\n"); - - _printf_(" parameters\n"); - if (parameters) parameters->Echo(); - else _printf_("parameters = NULL\n"); - - _printf_(" inputs\n"); - if (inputs) inputs->Echo(); - else _printf_("inputs=NULL\n"); } /*}}}*/ @@ -922,14 +809,4 @@ } /*}}}*/ -/*FUNCTION Tria::GetNumberOfNodesPressure;{{{*/ -int Tria::GetNumberOfNodesPressure(void){ - return this->NumberofNodesPressure(); -} -/*}}}*/ -/*FUNCTION Tria::GetNumberOfNodesVelocity;{{{*/ -int Tria::GetNumberOfNodesVelocity(void){ - return this->NumberofNodesVelocity(); -} -/*}}}*/ /*FUNCTION Tria::GetNumberOfVertices;{{{*/ int Tria::GetNumberOfVertices(void){ @@ -1008,18 +885,4 @@ return y; -} -/*}}}*/ -/*FUNCTION Tria::Id {{{*/ -int Tria::Id(){ - - return id; - -} -/*}}}*/ -/*FUNCTION Tria::Sid {{{*/ -int Tria::Sid(){ - - return sid; - } /*}}}*/ @@ -1351,17 +1214,4 @@ } /*}}}*/ -/*FUNCTION Tria::HasNodeOnBed {{{*/ -bool Tria::HasNodeOnBed(){ - - IssmDouble values[NUMVERTICES]; - IssmDouble sum; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&values[0],MeshVertexonbedEnum); - sum = values[0]+values[1]+values[2]; - - return sum>0.; -} -/*}}}*/ /*FUNCTION Tria::HasEdgeOnSurface {{{*/ bool Tria::HasEdgeOnSurface(){ @@ -1384,17 +1234,4 @@ return false; } -} -/*}}}*/ -/*FUNCTION Tria::HasNodeOnSurface {{{*/ -bool Tria::HasNodeOnSurface(){ - - IssmDouble values[NUMVERTICES]; - IssmDouble sum; - - /*Retrieve all inputs and parameters*/ - GetInputListOnVertices(&values[0],MeshVertexonbedEnum); - sum = values[0]+values[1]+values[2]; - - return sum>0.; } /*}}}*/ @@ -1540,18 +1377,4 @@ this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); return new GaussTria(indices[0],indices[1],order); -} -/*}}}*/ -/*FUNCTION Tria::IsIceInElement {{{*/ -bool Tria::IsIceInElement(){ - - /*Get levelset*/ - IssmDouble ls[NUMVERTICES]; - GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum); - - /*If the level set on one of the nodes is <0, ice is present in this element*/ - if(ls[0]<0. || ls[1]<0. || ls[2]<0.) - return true; - else - return false; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17516) @@ -32,7 +32,4 @@ public: - int id; - int sid; - /*Tria constructors, destructors {{{*/ Tria(){}; @@ -41,7 +38,4 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); int ObjectEnum(); Object *copy(); @@ -63,9 +57,6 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); - void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void Delta18oParameterization(void); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); - void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");}; - void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; int FiniteElement(void); Element* GetUpperElement(void){_error_("not implemented yet");}; @@ -77,14 +68,9 @@ int GetNodeIndex(Node* node); int GetNumberOfNodes(void); - int GetNumberOfNodesPressure(void); - int GetNumberOfNodesVelocity(void); int GetNumberOfVertices(void); - int Sid(); bool IsOnBed(); bool IsOnSurface(); bool HasEdgeOnBed(); - bool HasNodeOnBed(); bool HasEdgeOnSurface(); - bool HasNodeOnSurface(); void EdgeOnSurfaceIndices(int* pindex1,int* pindex); void EdgeOnBedIndices(int* pindex1,int* pindex); @@ -94,5 +80,4 @@ int NumberofNodesVelocity(void); int NumberofNodesPressure(void); - bool IsIceInElement(); void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); void GetVectorFromInputs(Vector* vector, int name_enum); @@ -111,5 +96,4 @@ Element* SpawnTopElement(void); int VelocityInterpolation(); - IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");}; int PressureInterpolation(); IssmDouble SurfaceArea(void); @@ -192,6 +176,4 @@ void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum); void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); - IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; - IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; IssmDouble GetArea(void); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17516) @@ -242,7 +242,4 @@ case HydrologyDCInefficientAnalysisEnum: Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax); - break; - case DamageEvolutionAnalysisEnum: - Ke=PenaltyCreateKMatrixDamageEvolution(kmax); break; default: @@ -276,7 +273,4 @@ case HydrologyDCInefficientAnalysisEnum: pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax); - break; - case DamageEvolutionAnalysisEnum: - pe=PenaltyCreatePVectorDamageEvolution(kmax); break; default: @@ -417,9 +411,4 @@ return; } - else if (analysis_type==DamageEvolutionAnalysisEnum){ - ConstraintActivateDamageEvolution(punstable); - return; - } - else{ _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet"); @@ -710,102 +699,4 @@ } /*}}}*/ -/*FUNCTION Pengrid::ConstraintActivateDamageEvolution {{{*/ -void Pengrid::ConstraintActivateDamageEvolution(int* punstable){ - - // The penalty is stable if it doesn't change during to successive iterations. - IssmDouble max_damage; - IssmDouble damage; - int new_active; - int unstable=0; - int penalty_lock; - - /*check that pengrid is not a clone (penalty to be added only once)*/ - if (node->IsClone()){ - unstable=0; - *punstable=unstable; - return; - } - - //First recover damage using the element: */ - element->GetInputValue(&damage,node,DamageDEnum); - - //Recover our data: - parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum); - parameters->FindParam(&max_damage,DamageMaxDamageEnum); - - //Figure out if damage>max_damage, in which case, this penalty needs to be activated. - //Would need to do the same for damage<0 if penalties are used. For now, ConstraintStatex - //is not called in solutionsequence_damage_nonlinear, so no penalties are applied. - - if (damage>max_damage){ - new_active=1; - } - else{ - new_active=0; - } - - //Figure out stability of this penalty - if (active==new_active){ - unstable=0; - } - else{ - unstable=1; - if(penalty_lock)zigzag_counter++; - } - - /*If penalty keeps zigzagging more than penalty_lock times: */ - if(penalty_lock){ - if(zigzag_counter>penalty_lock){ - unstable=0; - active=1; - } - } - - //Set penalty flag - active=new_active; - - //*Assign output pointers:*/ - *punstable=unstable; -} -/*}}}*/ -/*FUNCTION Pengrid::PenaltyCreateKMatrixDamageEvolution {{{*/ -ElementMatrix* Pengrid::PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax){ - - IssmDouble penalty_factor; - - /*Initialize Element matrix and return if necessary*/ - if(!this->active) return NULL; - ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters); - - /*recover parameters: */ - parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum); - - Ke->values[0]=kmax*pow(10.,penalty_factor); - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ -/*FUNCTION Pengrid::PenaltyCreatePVectorDamageEvolution {{{*/ -ElementVector* Pengrid::PenaltyCreatePVectorDamageEvolution(IssmDouble kmax){ - - IssmDouble penalty_factor; - IssmDouble max_damage; - - /*Initialize Element matrix and return if necessary*/ - if(!this->active) return NULL; - ElementVector* pe=new ElementVector(&node,1,this->parameters); - - //Recover our data: - parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum); - parameters->FindParam(&max_damage,DamageMaxDamageEnum); - - //right hand side penalizes to max_damage - pe->values[0]=kmax*pow(10.,penalty_factor)*max_damage; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ /*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/ ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){ Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17515) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17516) @@ -87,7 +87,4 @@ ElementVector* PenaltyCreatePVectorMelting(IssmDouble kmax); void ConstraintActivateThermal(int* punstable); - ElementMatrix* PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax); - ElementVector* PenaltyCreatePVectorDamageEvolution(IssmDouble kmax); - void ConstraintActivateDamageEvolution(int* punstable); ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax); ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax); Index: sm/trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp (revision 17515) +++ (revision ) @@ -1,80 +1,0 @@ -/* - * \brief: solutionsequence_damage_nonlinear.cpp: core of the damage solution - */ - -#include "../toolkits/toolkits.h" -#include "../classes/classes.h" -#include "../shared/shared.h" -#include "../modules/modules.h" - -void solutionsequence_damage_nonlinear(FemModel* femmodel){ - - /*solution : */ - Vector* Dg=NULL; - Vector* Df=NULL; - Vector* Df_old=NULL; - Vector* ys=NULL; - - /*intermediary: */ - Matrix* Kff=NULL; - Matrix* Kfs=NULL; - Vector* pf=NULL; - Vector* df=NULL; - - bool converged; - int constraints_converged; - int num_unstable_constraints; - int count; - int damage_penalty_threshold; - int damage_maxiter; - - /*parameters:*/ - int configuration_type; - - /*Recover parameters: */ - femmodel->parameters->FindParam(&damage_penalty_threshold,DamagePenaltyThresholdEnum); - femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); - femmodel->parameters->FindParam(&damage_maxiter,DamageMaxiterEnum); - - count=1; - converged=false; - - InputUpdateFromConstantx(femmodel,true,ResetPenaltiesEnum); - InputUpdateFromConstantx(femmodel,false,ConvergedEnum); - femmodel->UpdateConstraintsx(); - - for(;;){ - - delete Df_old; Df_old=Df; - SystemMatricesx(&Kff, &Kfs, &pf,&df, NULL,femmodel); - CreateNodalConstraintsx(&ys,femmodel->nodes,configuration_type); - Reduceloadx(pf, Kfs, ys); delete Kfs; - Solverx(&Df, Kff, pf,Df_old, df, femmodel->parameters); - delete Kff;delete pf;delete Dg; delete df; - Mergesolutionfromftogx(&Dg, Df,ys,femmodel->nodes,femmodel->parameters); delete ys; - InputUpdateFromSolutionx(femmodel,Dg); - - ConstraintsStatex(&constraints_converged,&num_unstable_constraints,femmodel); - - if (!converged){ - if(VerboseConvergence()) _printf0_(" #unstable constraints = " << num_unstable_constraints << "\n"); - if (num_unstable_constraints <= damage_penalty_threshold)converged=true; - if (count>=damage_maxiter){ - converged=true; - _printf0_(" maximum number of iterations (" << damage_maxiter << ") exceeded\n"); - } - } - count++; - - InputUpdateFromConstantx(femmodel,converged,ConvergedEnum); - - if(converged)break; - } - - InputUpdateFromSolutionx(femmodel,Dg); - - /*Free ressources: */ - delete Dg; - delete Df; - delete Df_old; -} Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17515) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17516) @@ -13,5 +13,4 @@ void solutionsequence_thermal_nonlinear(FemModel* femmodel); -void solutionsequence_damage_nonlinear(FemModel* femmodel); void solutionsequence_hydro_nonlinear(FemModel* femmodel); void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);