Index: ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17515) +++ ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp (revision 17516) @@ -213,8 +213,8 @@ IssmDouble *xyz_list = NULL; /*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; /*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/ @@ -319,8 +319,8 @@ } /*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*/ int* cs_list = xNew(vnumnodes+pnumnodes); @@ -938,8 +938,8 @@ IssmDouble FSreconditioning; /*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: ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17515) +++ ../trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17516) @@ -194,7 +194,7 @@ IssmDouble *vertex_pairing=NULL; 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;imy_vertices[reCast(vertex_pairing[2*i+1])-1]); /*Skip if one of the two is not on the bed*/ - if(!(reCast(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*/ penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); Index: ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 17515) +++ ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp (revision 17516) @@ -34,7 +34,7 @@ iomodel->FetchDataToInput(elements,SurfaceEnum); 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); iomodel->FetchDataToInput(elements,MeshElementonbedEnum); @@ -68,7 +68,7 @@ IssmDouble *vertex_pairing=NULL; 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;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ @@ -77,7 +77,9 @@ _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); /*Skip if one of the two is not on the bed*/ - if(!(reCast(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*/ penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17515) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17516) @@ -2753,8 +2753,8 @@ } /*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*/ ElementVector* de = element->NewElementVector(FSvelocityEnum); @@ -2783,8 +2783,8 @@ IssmDouble *xyz_list = NULL; /*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; /*Prepare coordinate system list*/ @@ -2885,8 +2885,8 @@ } /*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; @@ -2969,8 +2969,8 @@ } /*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; /*Initialize Element matrix and vectors*/ @@ -3064,8 +3064,8 @@ } /*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; /*Initialize Element matrix and vectors*/ @@ -3125,8 +3125,8 @@ } /*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*/ int* cs_list = xNew(vnumnodes+pnumnodes); @@ -3208,8 +3208,8 @@ } /*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*/ int* cs_list = xNew(vnumnodes+pnumnodes); @@ -3287,8 +3287,8 @@ } /*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*/ int* cs_list = xNew(vnumnodes+pnumnodes); @@ -3364,8 +3364,8 @@ } /*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(); /*Prepare coordinate system list*/ @@ -3658,8 +3658,8 @@ */ /*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; @@ -3785,8 +3785,8 @@ } /*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; @@ -4320,8 +4320,8 @@ element->GetInputValue(&approximation,ApproximationEnum); 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; /*Prepare node list*/ @@ -4438,8 +4438,8 @@ Element* pentabase=element->GetBasalElement(); 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; /*Prepare node list*/ @@ -4607,8 +4607,8 @@ 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(); int numnodes = vnumnodes+pnumnodes; /*Prepare coordinate system list*/ @@ -4686,8 +4686,8 @@ /*Initialize Element vector and return if necessary*/ 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*/ int* cs_list = xNew(vnumnodes+pnumnodes); @@ -4770,8 +4770,8 @@ if(!element->IsOnBed() || element->IsFloating()) return NULL; 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*/ int* cs_list = xNew(vnumnodes+pnumnodes); @@ -4846,8 +4846,8 @@ /*Initialize Element vector and return if necessary*/ 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*/ int* cs_list = xNew(vnumnodes+pnumnodes); Index: ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17515) +++ ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 17516) @@ -72,7 +72,7 @@ }/*}}}*/ 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); }/*}}}*/ @@ -86,21 +86,8 @@ }/*}}}*/ void DamageEvolutionAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ - /*create penalties for nodes: no node can have a damage > 1*/ - iomodel->FetchData(1,DamageSpcdamageEnum); - CreateSingleNodeToElementConnectivity(iomodel); + /*Nothing for now*/ - 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); - }/*}}}*/ /*Finite Element Analysis*/ Index: ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 17515) +++ ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp (revision 17516) @@ -59,7 +59,7 @@ IssmDouble *vertex_pairing=NULL; 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;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){ @@ -68,7 +68,9 @@ _assert_(iomodel->my_vertices[reCast(vertex_pairing[2*i+1])-1]); /*Skip if one of the two is not on the bed*/ - if(!(reCast(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*/ penpair_ids[0]=iomodel->nodecounter+reCast(vertex_pairing[2*i+0]); Index: ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp (revision 17515) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequence_damage_nonlinear.cpp (revision 17516) @@ -1,80 +0,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: ../trunk-jpl/src/c/solutionsequences/solutionsequences.h =================================================================== --- ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17515) +++ ../trunk-jpl/src/c/solutionsequences/solutionsequences.h (revision 17516) @@ -12,7 +12,6 @@ #include "../shared/Numerics/types.h" 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); void solutionsequence_newton(FemModel* femmodel); Index: ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17515) +++ ../trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 17516) @@ -242,9 +242,6 @@ case HydrologyDCInefficientAnalysisEnum: Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax); break; - case DamageEvolutionAnalysisEnum: - Ke=PenaltyCreateKMatrixDamageEvolution(kmax); - break; default: _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); } @@ -276,9 +273,6 @@ case HydrologyDCInefficientAnalysisEnum: pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax); break; - case DamageEvolutionAnalysisEnum: - pe=PenaltyCreatePVectorDamageEvolution(kmax); - break; default: _error_("analysis " << analysis_type << " (" << EnumToStringx(analysis_type) << ") not supported yet"); } @@ -416,11 +410,6 @@ ConstraintActivateHydrologyDCInefficient(punstable); return; } - else if (analysis_type==DamageEvolutionAnalysisEnum){ - ConstraintActivateDamageEvolution(punstable); - return; - } - else{ _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet"); } @@ -709,104 +698,6 @@ return pe; } /*}}}*/ -/*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: ../trunk-jpl/src/c/classes/Loads/Pengrid.h =================================================================== --- ../trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17515) +++ ../trunk-jpl/src/c/classes/Loads/Pengrid.h (revision 17516) @@ -86,9 +86,6 @@ ElementVector* PenaltyCreatePVectorThermal(IssmDouble kmax); 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); void ConstraintActivateHydrologyDCInefficient(int* punstable); Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17516) @@ -37,6 +37,8 @@ class Element: public Object,public Update{ public: + int id; + int sid; Inputs *inputs; Node **nodes; Vertex **vertices; @@ -54,8 +56,14 @@ /* bool AllActive(void); */ /* 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); void FindParam(IssmDouble* pvalue,int paramenum); @@ -80,19 +88,26 @@ void GetVerticesCoordinates(IssmDouble** xyz_list); 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); void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); + bool IsIceInElement(); bool IsInput(int name); bool IsFloating(); ElementVector* NewElementVector(int approximation_enum=NoneApproximationEnum); 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); void StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); @@ -134,12 +149,8 @@ virtual IssmDouble CharacteristicLength(void)=0; 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; virtual void NodalFunctions(IssmDouble* basis,Gauss* gauss)=0; @@ -153,7 +164,6 @@ virtual void NormalSection(IssmDouble* normal,IssmDouble* xyz_list)=0; 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; virtual Element* GetLowerElement(void)=0; @@ -168,16 +178,11 @@ virtual void GetSolutionFromInputsOneDof(Vector* solution,int solutionenum)=0; 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; virtual void GetInputValue(IssmDouble* pvalue,Node* node,int enumtype)=0; Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17516) @@ -94,64 +94,6 @@ /*}}}*/ /*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){ @@ -291,35 +233,6 @@ } /*}}}*/ -/*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){ @@ -415,32 +328,6 @@ *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"); -} -/*}}}*/ /*FUNCTION Tria::FiniteElement{{{*/ int Tria::FiniteElement(void){ return this->element_type; @@ -921,16 +808,6 @@ return this->NumberofNodes(); } /*}}}*/ -/*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){ return NUMVERTICES; @@ -1009,20 +886,6 @@ return y; } /*}}}*/ -/*FUNCTION Tria::Id {{{*/ -int Tria::Id(){ - - return id; - -} -/*}}}*/ -/*FUNCTION Tria::Sid {{{*/ -int Tria::Sid(){ - - return sid; - -} -/*}}}*/ /*FUNCTION Tria::InputDepthAverageAtBase {{{*/ void Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){ @@ -1350,19 +1213,6 @@ } } /*}}}*/ -/*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(){ @@ -1385,19 +1235,6 @@ } } /*}}}*/ -/*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.; -} -/*}}}*/ /*FUNCTION Tria::EdgeOnBedIndices{{{*/ void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){ @@ -1541,20 +1378,6 @@ 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; -} -/*}}}*/ /*FUNCTION Tria::NodalFunctions{{{*/ void Tria::NodalFunctions(IssmDouble* basis, Gauss* gauss){ Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17516) @@ -31,18 +31,12 @@ public: - int id; - int sid; - /*Tria constructors, destructors {{{*/ Tria(){}; Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels); ~Tria(); /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); int ObjectEnum(); Object *copy(); /*}}}*/ @@ -62,11 +56,8 @@ void ComputeSurfaceNormalVelocity(); 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");}; Element* GetLowerElement(void){_error_("not implemented yet");}; @@ -76,16 +67,11 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); 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); int EdgeOnBedIndex(); @@ -93,7 +79,6 @@ bool IsNodeOnShelfFromFlags(IssmDouble* flags); int NumberofNodesVelocity(void); int NumberofNodesPressure(void); - bool IsIceInElement(); void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); void GetVectorFromInputs(Vector* vector, int name_enum); void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); @@ -110,7 +95,6 @@ Element* SpawnBasalElement(void); Element* SpawnTopElement(void); int VelocityInterpolation(); - IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");}; int PressureInterpolation(); IssmDouble SurfaceArea(void); void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement); @@ -191,8 +175,6 @@ /*Tria specific routines:{{{*/ 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); int GetElementType(void); Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17516) @@ -370,26 +370,6 @@ this->inputs->Configure(parameters); } /*}}}*/ -/*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){ /*Are we on the base? If not, return*/ @@ -464,32 +444,6 @@ 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); -} -/*}}}*/ /*FUNCTION Penta::FiniteElement{{{*/ int Penta::FiniteElement(void){ return this->element_type; @@ -847,16 +801,6 @@ return this->NumberofNodes(); } /*}}}*/ -/*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){ return NUMVERTICES; @@ -1161,18 +1105,6 @@ *pxyz_zero= xyz_zero; } /*}}}*/ -/*FUNCTION Penta::Sid {{{*/ -int Penta::Sid(){ - - return sid; - -} -/*}}}*/ -/*FUNCTION Penta::Id {{{*/ -int Penta::Id(void){ - return id; -} -/*}}}*/ /*FUNCTION Penta::InputDepthAverageAtBase{{{*/ void Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){ @@ -1642,20 +1574,6 @@ } /*}}}*/ -/*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){ /*Return the minimum lenght of the nine egdes of the penta*/ @@ -1812,12 +1730,11 @@ /*FUNCTION Penta::NormalBase {{{*/ 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]; } @@ -1825,7 +1742,7 @@ normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; 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*/ bed_normal[0]=-normal[0]/normal_norm; @@ -1918,12 +1835,6 @@ delete gauss; } /*}}}*/ -/*FUNCTION Penta::PureIceEnthalpy{{{*/ -IssmDouble Penta::PureIceEnthalpy(IssmDouble pressure){ - - return this->matpar->PureIceEnthalpy(pressure); -} -/*}}}*/ /*FUNCTION Penta::ReduceMatrices{{{*/ void Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){ @@ -2071,77 +1982,6 @@ 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; -} -/*}}}*/ /*FUNCTION Penta::SpawnTria {{{*/ Tria* Penta::SpawnTria(int index1,int index2,int index3){ Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17516) @@ -31,9 +31,6 @@ public: - int id; - int sid; - Penta **verticalneighbors; // 2 neighbors: first one under, second one above /*Penta constructors and destructor: {{{*/ @@ -43,10 +40,7 @@ /*}}}*/ /*Object virtual functions definitions: {{{*/ Object *copy(); - void DeepEcho(); - void Echo(); int ObjectEnum(); - int Id(); /*}}}*/ /*Update virtual functions definitions: {{{*/ void InputUpdateFromVector(IssmDouble* vector, int name, int type); @@ -66,10 +60,7 @@ void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); 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); Penta* GetSurfacePenta(void); @@ -82,8 +73,6 @@ IssmDouble GetGroundedPortion(IssmDouble* xyz_list); int GetNodeIndex(Node* node); int GetNumberOfNodes(void); - int GetNumberOfNodesPressure(void); - int GetNumberOfNodesVelocity(void); int GetNumberOfVertices(void); void GetSolutionFromInputsOneDof(Vector* solution,int enum_type); IssmDouble GetXcoord(Gauss* gauss); @@ -93,14 +82,12 @@ void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); - int Sid(); void InputDepthAverageAtBase(int enum_type,int average_enum_type); void InputDuplicate(int original_enum,int new_enum); void InputScale(int enum_type,IssmDouble scale_factor); int NumberofNodesVelocity(void); int NumberofNodesPressure(void); int VelocityInterpolation(); - IssmDouble PureIceEnthalpy(IssmDouble pressure); int PressureInterpolation(); bool IsZeroLevelset(int levelset_enum); bool IsIcefront(void); @@ -188,8 +175,6 @@ void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); 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); void GetVertexPidList(int* doflist); @@ -207,7 +192,6 @@ void JacobianDeterminantSurface(IssmDouble* pJdet, IssmDouble* xyz_list,Gauss* gauss); 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); Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");}; Index: ../trunk-jpl/src/c/classes/Elements/Seg.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Seg.cpp (revision 17516) @@ -64,64 +64,11 @@ return sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)); } /*}}}*/ -/*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){ return this->element_type; } /*}}}*/ -/*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){ @@ -129,14 +76,7 @@ } /*}}}*/ -/*FUNCTION Seg::Id {{{*/ -int Seg::Id(){ - return id; - -} -/*}}}*/ - void Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/ /* Intermediaries */ @@ -185,20 +125,6 @@ *pxyz_list = xyz_list; }/*}}}*/ -/*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){/*{{{*/ bool isicefront; Index: ../trunk-jpl/src/c/classes/Elements/Tetra.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Tetra.cpp (revision 17516) @@ -51,80 +51,16 @@ } /*}}}*/ -/*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){ return TetraEnum; -} -/*}}}*/ -/*FUNCTION Tetra::Id {{{*/ -int Tetra::Id(){ - return id; -} -/*}}}*/ +}/*}}}*/ /*FUNCTION Tetra::AddInput{{{*/ void Tetra::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){ @@ -239,16 +175,6 @@ return NUMVERTICES; } /*}}}*/ -/*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){ @@ -340,19 +266,6 @@ } } /*}}}*/ -/*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.; -} -/*}}}*/ /*FUNCTION Tetra::InputUpdateFromIoModel {{{*/ void Tetra::InputUpdateFromIoModel(int index,IoModel* iomodel){ @@ -484,21 +397,6 @@ xDelete(doflist); } /*}}}*/ -/*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){ int size = nodes[4]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); @@ -680,7 +624,8 @@ 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){ int size = nodes[4]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum); @@ -770,71 +715,6 @@ } /*}}}*/ -/*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; - -} -/*}}}*/ /*FUNCTION Tetra::SpawnBasalElement{{{*/ Element* Tetra::SpawnBasalElement(void){ Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17516) @@ -29,18 +29,12 @@ public: - int id; - int sid; - /*Seg constructors, destructors {{{*/ Seg(){}; Seg(int seg_id,int seg_sid,int i, IoModel* iomodel,int nummodels); ~Seg(); /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); int ObjectEnum(); Object *copy(); /*}}}*/ @@ -63,13 +57,8 @@ 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 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");}; Element* GetLowerElement(void){_error_("not implemented yet");}; @@ -77,13 +66,10 @@ Element* GetBasalElement(void){_error_("not implemented yet");}; 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");}; bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; @@ -101,7 +87,6 @@ void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; 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");}; void NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; @@ -110,7 +95,6 @@ Element* SpawnBasalElement(void){_error_("not implemented yet");}; 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");}; void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 17516) @@ -29,18 +29,12 @@ public: - int id; - int sid; - /*Tetra constructors, destructors {{{*/ Tetra(){}; Tetra(int seg_id,int seg_sid,int i, IoModel* iomodel,int nummodels); ~Tetra(); /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); int ObjectEnum(); Object *copy(); /*}}}*/ @@ -63,13 +57,8 @@ void ComputeStressTensor(){_error_("not implemented yet");}; 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); void FaceOnSurfaceIndices(int* pindex1,int* pindex2,int* pindex3); @@ -80,15 +69,11 @@ Element* GetBasalElement(void){_error_("not implemented yet");}; 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); void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); bool HasFaceOnBed(); bool HasFaceOnSurface(); - bool HasNodeOnBed(); - int Sid(); bool IsOnBed(); bool IsOnSurface(); bool IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");}; @@ -106,17 +91,15 @@ void NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; 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); Element* SpawnBasalElement(void); Element* SpawnTopElement(void); 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");}; void ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; @@ -227,4 +210,4 @@ int UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){_error_("not implemented yet");}; /*}}}*/ }; -#endif /* _SEG_H */ +#endif /* _TETRA_H_*/ Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17515) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 17516) @@ -16,6 +16,8 @@ /*Constructors/destructor/copy*/ Element::Element(){/*{{{*/ + this->id = -1; + this->sid = -1; this->inputs = NULL; this->nodes = NULL; this->vertices = NULL; @@ -129,6 +131,71 @@ void Element::DeleteMaterials(void){/*{{{*/ 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*/ @@ -161,6 +228,18 @@ delete gauss; 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); }/*}}}*/ @@ -199,7 +278,7 @@ void Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/ /*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: */ int numberofdofs=0; @@ -222,8 +301,8 @@ void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/ /*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: */ int numberofdofs=0; @@ -387,7 +466,7 @@ _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"); @@ -469,6 +548,18 @@ 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; + +} +/*}}}*/ void Element::InputChangeName(int original_enum,int new_enum){/*{{{*/ this->inputs->ChangeEnum(original_enum,new_enum); } @@ -586,6 +677,10 @@ return shelf; }/*}}}*/ +bool Element::IsIceInElement(){/*{{{*/ + return (this->inputs->Min(MaskIceLevelsetEnum)<0.); +} +/*}}}*/ bool Element::IsInput(int name){/*{{{*/ if ( name==ThicknessEnum || @@ -674,6 +769,9 @@ return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum); } /*}}}*/ +IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/ + return this->matpar->PureIceEnthalpy(pressure); +}/*}}}*/ void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/ Input* input=this->inputs->GetInput(output_enum); @@ -772,6 +870,82 @@ input->ResultToPatch(values,nodesperelement,this->Sid()); } /*}}}*/ +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); return this->matpar->TMeltingPoint(pressure);