![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 17516
- Timestamp:
- 03/21/14 21:07:11 (11 years ago)
- Location:
- issm/trunk-jpl/src/c
- Files:
-
- 1 deleted
- 19 edited
-
analyses/AdjointHorizAnalysis.cpp (modified) (3 diffs)
-
analyses/DamageEvolutionAnalysis.cpp (modified) (2 diffs)
-
analyses/FreeSurfaceBaseAnalysis.cpp (modified) (2 diffs)
-
analyses/FreeSurfaceTopAnalysis.cpp (modified) (3 diffs)
-
analyses/MasstransportAnalysis.cpp (modified) (2 diffs)
-
analyses/StressbalanceAnalysis.cpp (modified) (17 diffs)
-
classes/Elements/Element.cpp (modified) (10 diffs)
-
classes/Elements/Element.h (modified) (8 diffs)
-
classes/Elements/Penta.cpp (modified) (9 diffs)
-
classes/Elements/Penta.h (modified) (8 diffs)
-
classes/Elements/Seg.cpp (modified) (3 diffs)
-
classes/Elements/Seg.h (modified) (6 diffs)
-
classes/Elements/Tetra.cpp (modified) (9 diffs)
-
classes/Elements/Tetra.h (modified) (8 diffs)
-
classes/Elements/Tria.cpp (modified) (8 diffs)
-
classes/Elements/Tria.h (modified) (7 diffs)
-
classes/Loads/Pengrid.cpp (modified) (4 diffs)
-
classes/Loads/Pengrid.h (modified) (1 diff)
-
solutionsequences/solutionsequence_damage_nonlinear.cpp (deleted)
-
solutionsequences/solutionsequences.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp
r17212 r17516 214 214 215 215 /*Fetch number of nodes and dof for this finite element*/ 216 int vnumnodes = element-> GetNumberOfNodesVelocity();217 int pnumnodes = element-> GetNumberOfNodesPressure();216 int vnumnodes = element->NumberofNodesVelocity(); 217 int pnumnodes = element->NumberofNodesPressure(); 218 218 int numdof = vnumnodes*3 + pnumnodes; 219 219 … … 320 320 321 321 /*Fetch number of nodes and dof for this finite element*/ 322 int vnumnodes = element-> GetNumberOfNodesVelocity();323 int pnumnodes = element-> GetNumberOfNodesPressure();322 int vnumnodes = element->NumberofNodesVelocity(); 323 int pnumnodes = element->NumberofNodesPressure(); 324 324 325 325 /*Prepare coordinate system list*/ … … 939 939 940 940 /*Fetch number of nodes and dof for this finite element*/ 941 int vnumnodes = element-> GetNumberOfNodesVelocity();942 int pnumnodes = element-> GetNumberOfNodesPressure();941 int vnumnodes = element->NumberofNodesVelocity(); 942 int pnumnodes = element->NumberofNodesPressure(); 943 943 int vnumdof = vnumnodes*3; 944 944 int pnumdof = pnumnodes*1; -
issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp
r17396 r17516 73 73 void DamageEvolutionAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/ 74 74 75 i omodel->FetchData(1,MeshVertexonbedEnum);75 if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(1,MeshVertexonbedEnum); 76 76 ::CreateNodes(nodes,iomodel,DamageEvolutionAnalysisEnum,P1Enum); 77 77 iomodel->DeleteData(1,MeshVertexonbedEnum); … … 87 87 void DamageEvolutionAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ 88 88 89 /*create penalties for nodes: no node can have a damage > 1*/ 90 iomodel->FetchData(1,DamageSpcdamageEnum); 91 CreateSingleNodeToElementConnectivity(iomodel); 92 93 for(int i=0;i<iomodel->numberofvertices;i++){ 94 95 /*keep only this partition's nodes:*/ 96 if(iomodel->my_vertices[i]){ 97 if (xIsNan<IssmDouble>(iomodel->Data(DamageSpcdamageEnum)[i])){ //No penalty applied on spc nodes! 98 loads->AddObject(new Pengrid(iomodel->loadcounter+i+1,i,iomodel,DamageEvolutionAnalysisEnum)); 99 } 100 } 101 } 102 iomodel->DeleteData(1,DamageSpcdamageEnum); 89 /*Nothing for now*/ 103 90 104 91 }/*}}}*/ -
issm/trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
r17212 r17516 60 60 IssmDouble *nodeonbed=NULL; 61 61 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); 62 i omodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);62 if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); 63 63 for(int i=0;i<numvertex_pairing;i++){ 64 64 … … 69 69 70 70 /*Skip if one of the two is not on the bed*/ 71 if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue; 71 if(iomodel->meshtype!=Mesh2DhorizontalEnum){ 72 if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue; 73 } 72 74 73 75 /*Get node ids*/ -
issm/trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp
r17212 r17516 35 35 iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum); 36 36 iomodel->FetchDataToInput(elements,VxEnum); 37 i omodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);37 if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum); 38 38 if(iomodel->meshtype==Mesh3DEnum){ 39 39 iomodel->FetchDataToInput(elements,VzEnum); … … 69 69 IssmDouble *nodeonsurface=NULL; 70 70 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); 71 i omodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);71 if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum); 72 72 for(int i=0;i<numvertex_pairing;i++){ 73 73 … … 78 78 79 79 /*Skip if one of the two is not on the bed*/ 80 if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue; 80 if(iomodel->meshtype!=Mesh2DhorizontalEnum){ 81 if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue; 82 } 81 83 82 84 /*Get node ids*/ -
issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp
r17463 r17516 195 195 IssmDouble *nodeonbed=NULL; 196 196 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum); 197 i omodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum);197 if(iomodel->meshtype!=Mesh2DhorizontalEnum) iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); 198 198 199 199 for(int i=0;i<numvertex_pairing;i++){ … … 205 205 206 206 /*Skip if one of the two is not on the bed*/ 207 if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue; 207 if(iomodel->meshtype!=Mesh2DhorizontalEnum){ 208 if(!(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonbed[reCast<int>(vertex_pairing[2*i+1])-1]))) continue; 209 } 208 210 209 211 /*Get node ids*/ -
issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
r17493 r17516 2754 2754 2755 2755 /*Fetch number of nodes and dof for this finite element*/ 2756 int vnumnodes = element-> GetNumberOfNodesVelocity();2757 int pnumnodes = element-> GetNumberOfNodesPressure();2756 int vnumnodes = element->NumberofNodesVelocity(); 2757 int pnumnodes = element->NumberofNodesPressure(); 2758 2758 2759 2759 /*Initialize output vector*/ … … 2784 2784 2785 2785 /*Fetch number of nodes and dof for this finite element*/ 2786 int vnumnodes = element-> GetNumberOfNodesVelocity();2787 int pnumnodes = element-> GetNumberOfNodesPressure();2786 int vnumnodes = element->NumberofNodesVelocity(); 2787 int pnumnodes = element->NumberofNodesPressure(); 2788 2788 int numdof = vnumnodes*NDOF3 + pnumnodes*NDOF1; 2789 2789 … … 2886 2886 2887 2887 /*Fetch number of nodes and dof for this finite element*/ 2888 int vnumnodes = element-> GetNumberOfNodesVelocity();2889 int pnumnodes = element-> GetNumberOfNodesPressure();2888 int vnumnodes = element->NumberofNodesVelocity(); 2889 int pnumnodes = element->NumberofNodesPressure(); 2890 2890 int numdof = vnumnodes*dim + pnumnodes; 2891 2891 int bsize = epssize + 2; … … 2970 2970 2971 2971 /*Fetch number of nodes and dof for this finite element*/ 2972 int vnumnodes = element-> GetNumberOfNodesVelocity();2973 int pnumnodes = element-> GetNumberOfNodesPressure();2972 int vnumnodes = element->NumberofNodesVelocity(); 2973 int pnumnodes = element->NumberofNodesPressure(); 2974 2974 int numdof = vnumnodes*dim + pnumnodes; 2975 2975 … … 3065 3065 3066 3066 /*Fetch number of nodes and dof for this finite element*/ 3067 int vnumnodes = element-> GetNumberOfNodesVelocity();3068 int pnumnodes = element-> GetNumberOfNodesPressure();3067 int vnumnodes = element->NumberofNodesVelocity(); 3068 int pnumnodes = element->NumberofNodesPressure(); 3069 3069 int numdof = vnumnodes*dim + pnumnodes; 3070 3070 … … 3126 3126 3127 3127 /*Fetch number of nodes and dof for this finite element*/ 3128 int vnumnodes = element-> GetNumberOfNodesVelocity();3129 int pnumnodes = element-> GetNumberOfNodesPressure();3128 int vnumnodes = element->NumberofNodesVelocity(); 3129 int pnumnodes = element->NumberofNodesPressure(); 3130 3130 3131 3131 /*Prepare coordinate system list*/ … … 3209 3209 3210 3210 /*Fetch number of nodes and dof for this finite element*/ 3211 int vnumnodes = element-> GetNumberOfNodesVelocity();3212 int pnumnodes = element-> GetNumberOfNodesPressure();3211 int vnumnodes = element->NumberofNodesVelocity(); 3212 int pnumnodes = element->NumberofNodesPressure(); 3213 3213 3214 3214 /*Prepare coordinate system list*/ … … 3288 3288 3289 3289 /*Fetch number of nodes and dof for this finite element*/ 3290 int vnumnodes = element-> GetNumberOfNodesVelocity();3291 int pnumnodes = element-> GetNumberOfNodesPressure();3290 int vnumnodes = element->NumberofNodesVelocity(); 3291 int pnumnodes = element->NumberofNodesPressure(); 3292 3292 3293 3293 /*Prepare coordinate system list*/ … … 3365 3365 3366 3366 /*Fetch number of nodes and dof for this finite element*/ 3367 int vnumnodes = element-> GetNumberOfNodesVelocity();3368 int pnumnodes = element-> GetNumberOfNodesPressure();3367 int vnumnodes = element->NumberofNodesVelocity(); 3368 int pnumnodes = element->NumberofNodesPressure(); 3369 3369 int numvertices = element->GetNumberOfVertices(); 3370 3370 … … 3659 3659 3660 3660 /*Fetch number of nodes for this finite element*/ 3661 int pnumnodes = element-> GetNumberOfNodesPressure();3662 int vnumnodes = element-> GetNumberOfNodesVelocity();3661 int pnumnodes = element->NumberofNodesPressure(); 3662 int vnumnodes = element->NumberofNodesVelocity(); 3663 3663 int pnumdof = pnumnodes; 3664 3664 int vnumdof = vnumnodes*dim; … … 3786 3786 3787 3787 /*Fetch number of nodes and dof for this finite element*/ 3788 int vnumnodes = element-> GetNumberOfNodesVelocity();3789 int pnumnodes = element-> GetNumberOfNodesPressure();3788 int vnumnodes = element->NumberofNodesVelocity(); 3789 int pnumnodes = element->NumberofNodesPressure(); 3790 3790 int vnumdof = vnumnodes*dim; 3791 3791 int pnumdof = pnumnodes*1; … … 4321 4321 if(element->IsFloating() || !element->IsOnBed()) return NULL; 4322 4322 4323 int vnumnodes = element-> GetNumberOfNodesVelocity();4324 int pnumnodes = element-> GetNumberOfNodesPressure();4323 int vnumnodes = element->NumberofNodesVelocity(); 4324 int pnumnodes = element->NumberofNodesPressure(); 4325 4325 int numnodes = 2*vnumnodes-1+pnumnodes; 4326 4326 … … 4439 4439 Element* basaltria=pentabase->SpawnBasalElement(); 4440 4440 4441 int vnumnodes = element-> GetNumberOfNodesVelocity();4442 int pnumnodes = element-> GetNumberOfNodesPressure();4441 int vnumnodes = element->NumberofNodesVelocity(); 4442 int pnumnodes = element->NumberofNodesPressure(); 4443 4443 int numnodes = 2*vnumnodes-1+pnumnodes; 4444 4444 … … 4608 4608 if(approximation!=HOFSApproximationEnum) return NULL; 4609 4609 4610 int vnumnodes = element-> GetNumberOfNodesVelocity();4611 int pnumnodes = element-> GetNumberOfNodesPressure();4610 int vnumnodes = element->NumberofNodesVelocity(); 4611 int pnumnodes = element->NumberofNodesPressure(); 4612 4612 int numnodes = vnumnodes+pnumnodes; 4613 4613 … … 4687 4687 element->GetInputValue(&approximation,ApproximationEnum); 4688 4688 if(approximation!=HOFSApproximationEnum) return NULL; 4689 int vnumnodes = element-> GetNumberOfNodesVelocity();4690 int pnumnodes = element-> GetNumberOfNodesPressure();4689 int vnumnodes = element->NumberofNodesVelocity(); 4690 int pnumnodes = element->NumberofNodesPressure(); 4691 4691 4692 4692 /*Prepare coordinate system list*/ … … 4771 4771 element->GetInputValue(&approximation,ApproximationEnum); 4772 4772 if(approximation!=SSAFSApproximationEnum) return NULL; 4773 int vnumnodes = element-> GetNumberOfNodesVelocity();4774 int pnumnodes = element-> GetNumberOfNodesPressure();4773 int vnumnodes = element->NumberofNodesVelocity(); 4774 int pnumnodes = element->NumberofNodesPressure(); 4775 4775 4776 4776 /*Prepare coordinate system list*/ … … 4847 4847 element->GetInputValue(&approximation,ApproximationEnum); 4848 4848 if(approximation!=SSAFSApproximationEnum) return NULL; 4849 int vnumnodes = element-> GetNumberOfNodesVelocity();4850 int pnumnodes = element-> GetNumberOfNodesPressure();4849 int vnumnodes = element->NumberofNodesVelocity(); 4850 int pnumnodes = element->NumberofNodesPressure(); 4851 4851 4852 4852 /*Prepare coordinate system list*/ -
issm/trunk-jpl/src/c/classes/Elements/Element.cpp
r17472 r17516 17 17 /*Constructors/destructor/copy*/ 18 18 Element::Element(){/*{{{*/ 19 this->id = -1; 20 this->sid = -1; 19 21 this->inputs = NULL; 20 22 this->nodes = NULL; … … 130 132 delete this->material; 131 133 }/*}}}*/ 134 void Element::DeepEcho(void){/*{{{*/ 135 136 _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n"); 137 _printf_(" id : "<<this->id <<"\n"); 138 _printf_(" sid: "<<this->sid<<"\n"); 139 if(vertices){ 140 int numvertices = this->GetNumberOfVertices(); 141 for(int i=0;i<numvertices;i++) vertices[i]->Echo(); 142 } 143 else _printf_("vertices = NULL\n"); 144 145 if(nodes){ 146 int numnodes = this->GetNumberOfNodes(); 147 for(int i=0;i<numnodes;i++) nodes[i]->DeepEcho(); 148 } 149 else _printf_("nodes = NULL\n"); 150 151 if (material) material->DeepEcho(); 152 else _printf_("material = NULL\n"); 153 154 if (matpar) matpar->DeepEcho(); 155 else _printf_("matpar = NULL\n"); 156 157 _printf_(" parameters\n"); 158 if (parameters) parameters->DeepEcho(); 159 else _printf_("parameters = NULL\n"); 160 161 _printf_(" inputs\n"); 162 if (inputs) inputs->DeepEcho(); 163 else _printf_("inputs=NULL\n"); 164 165 return; 166 } 167 /*}}}*/ 168 void Element::Echo(void){/*{{{*/ 169 _printf_(EnumToStringx(this->ObjectEnum())<<" element:\n"); 170 _printf_(" id : "<<this->id <<"\n"); 171 _printf_(" sid: "<<this->sid<<"\n"); 172 if(vertices){ 173 int numvertices = this->GetNumberOfVertices(); 174 for(int i=0;i<numvertices;i++) vertices[i]->Echo(); 175 } 176 else _printf_("vertices = NULL\n"); 177 178 if(nodes){ 179 int numnodes = this->GetNumberOfNodes(); 180 for(int i=0;i<numnodes;i++) nodes[i]->Echo(); 181 } 182 else _printf_("nodes = NULL\n"); 183 184 if (material) material->Echo(); 185 else _printf_("material = NULL\n"); 186 187 if (matpar) matpar->Echo(); 188 else _printf_("matpar = NULL\n"); 189 190 _printf_(" parameters\n"); 191 if (parameters) parameters->Echo(); 192 else _printf_("parameters = NULL\n"); 193 194 _printf_(" inputs\n"); 195 if (inputs) inputs->Echo(); 196 else _printf_("inputs=NULL\n"); 197 } 198 /*}}}*/ 132 199 IssmDouble Element::Divergence(void){/*{{{*/ 133 200 /*Compute element divergence*/ … … 162 229 return divergence; 163 230 }/*}}}*/ 231 void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/ 232 matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure); 233 }/*}}}*/ 234 void Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ 235 matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure); 236 }/*}}}*/ 237 IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/ 238 return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); 239 }/*}}}*/ 240 IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/ 241 return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); 242 }/*}}}*/ 164 243 void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/ 165 244 this->parameters->FindParam(pvalue,paramenum); … … 200 279 201 280 /*Fetch number of nodes and dof for this finite element*/ 202 int numnodes = this-> GetNumberOfNodesVelocity();281 int numnodes = this->NumberofNodesVelocity(); 203 282 204 283 /*First, figure out size of doflist and create it: */ … … 223 302 224 303 /*Fetch number of nodes and dof for this finite element*/ 225 int vnumnodes = this-> GetNumberOfNodesVelocity();226 int pnumnodes = this-> GetNumberOfNodesPressure();304 int vnumnodes = this->NumberofNodesVelocity(); 305 int pnumnodes = this->NumberofNodesPressure(); 227 306 228 307 /*First, figure out size of doflist and create it: */ … … 388 467 _assert_(pvalue); 389 468 390 int numnodes = this-> GetNumberOfNodesVelocity();469 int numnodes = this->NumberofNodesVelocity(); 391 470 Input *input = this->GetInput(enumtype); 392 471 if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element"); … … 468 547 int numvertices = this->GetNumberOfVertices(); 469 548 for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity(); 549 } 550 /*}}}*/ 551 bool Element::HasNodeOnBed(){/*{{{*/ 552 return (this->inputs->Max(MeshVertexonbedEnum)>0.); 553 }/*}}}*/ 554 bool Element::HasNodeOnSurface(){/*{{{*/ 555 return (this->inputs->Max(MeshVertexonsurfaceEnum)>0.); 556 }/*}}}*/ 557 int Element::Id(){/*{{{*/ 558 559 return this->id; 560 470 561 } 471 562 /*}}}*/ … … 587 678 return shelf; 588 679 }/*}}}*/ 680 bool Element::IsIceInElement(){/*{{{*/ 681 return (this->inputs->Min(MaskIceLevelsetEnum)<0.); 682 } 683 /*}}}*/ 589 684 bool Element::IsInput(int name){/*{{{*/ 590 685 if ( … … 675 770 } 676 771 /*}}}*/ 772 IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/ 773 return this->matpar->PureIceEnthalpy(pressure); 774 }/*}}}*/ 677 775 void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/ 678 776 … … 773 871 774 872 } /*}}}*/ 873 void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/ 874 875 /*Intermediaries*/ 876 const int numnodes = this->GetNumberOfNodes(); 877 878 /*Output */ 879 int d_nz = 0; 880 int o_nz = 0; 881 882 /*Loop over all nodes*/ 883 for(int i=0;i<numnodes;i++){ 884 885 if(!flags[this->nodes[i]->Lid()]){ 886 887 /*flag current node so that no other element processes it*/ 888 flags[this->nodes[i]->Lid()]=true; 889 890 int counter=0; 891 while(flagsindices[counter]>=0) counter++; 892 flagsindices[counter]=this->nodes[i]->Lid(); 893 894 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/ 895 switch(set2_enum){ 896 case FsetEnum: 897 if(nodes[i]->indexing.fsize){ 898 if(this->nodes[i]->IsClone()) 899 o_nz += 1; 900 else 901 d_nz += 1; 902 } 903 break; 904 case GsetEnum: 905 if(nodes[i]->indexing.gsize){ 906 if(this->nodes[i]->IsClone()) 907 o_nz += 1; 908 else 909 d_nz += 1; 910 } 911 break; 912 case SsetEnum: 913 if(nodes[i]->indexing.ssize){ 914 if(this->nodes[i]->IsClone()) 915 o_nz += 1; 916 else 917 d_nz += 1; 918 } 919 break; 920 default: _error_("not supported"); 921 } 922 } 923 } 924 925 /*Special case: 2d/3d coupling, the node of this element might be connected 926 *to the basal element*/ 927 int analysis_type,approximation,numlayers; 928 parameters->FindParam(&analysis_type,AnalysisTypeEnum); 929 if(analysis_type==StressbalanceAnalysisEnum){ 930 inputs->GetInputValue(&approximation,ApproximationEnum); 931 if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){ 932 parameters->FindParam(&numlayers,MeshNumberoflayersEnum); 933 o_nz += numlayers*3; 934 d_nz += numlayers*3; 935 } 936 } 937 938 /*Assign output pointers: */ 939 *pd_nz=d_nz; 940 *po_nz=o_nz; 941 } 942 /*}}}*/ 943 int Element::Sid(){/*{{{*/ 944 945 return this->sid; 946 947 } 948 /*}}}*/ 775 949 IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/ 776 950 _assert_(matpar); -
issm/trunk-jpl/src/c/classes/Elements/Element.h
r17472 r17516 38 38 39 39 public: 40 int id; 41 int sid; 40 42 Inputs *inputs; 41 43 Node **nodes; … … 55 57 /* bool AnyActive(void); */ 56 58 void CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes,int numnodes,int* cs_array); 59 void Echo(); 60 void DeepEcho(); 57 61 void DeleteMaterials(void); 58 62 IssmDouble Divergence(void); 63 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure); 64 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); 65 IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); 66 IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); 59 67 void FindParam(bool* pvalue,int paramenum); 60 68 void FindParam(int* pvalue,int paramenum); … … 81 89 void GetVerticesSidList(int* sidlist); 82 90 void GetVerticesConnectivityList(int* connectivitylist); 91 bool HasNodeOnBed(); 92 bool HasNodeOnSurface(); 93 int Id(); 94 int Sid(); 83 95 void InputChangeName(int enum_type,int enum_type_old); 84 96 void InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code); … … 86 98 void InputUpdateFromConstant(int constant, int name); 87 99 void InputUpdateFromConstant(bool constant, int name); 100 bool IsIceInElement(); 88 101 bool IsInput(int name); 89 102 bool IsFloating(); … … 91 104 ElementMatrix* NewElementMatrix(int approximation_enum=NoneApproximationEnum); 92 105 ElementMatrix* NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum); 106 IssmDouble PureIceEnthalpy(IssmDouble pressure); 93 107 void ResultInterpolation(int* pinterpolation,int*nodesperelement,int output_enum); 94 108 void ResultToVector(Vector<IssmDouble>* vector,int output_enum); 95 109 void ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum); 110 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); 96 111 void StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input); 97 112 void StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input); … … 135 150 virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; 136 151 virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters)=0; 137 virtual void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0;138 152 virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; 139 virtual void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure)=0; 140 virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0; 141 virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0; 142 virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0; 153 143 154 virtual int FiniteElement(void)=0; 144 155 virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0; … … 154 165 virtual void NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0; 155 166 virtual void NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0; 156 virtual IssmDouble PureIceEnthalpy(IssmDouble pressure)=0;157 167 158 168 virtual Element* GetUpperElement(void)=0; … … 169 179 virtual int GetNodeIndex(Node* node)=0; 170 180 virtual int GetNumberOfNodes(void)=0; 171 virtual int GetNumberOfNodesVelocity(void)=0;172 virtual int GetNumberOfNodesPressure(void)=0;173 181 virtual int GetNumberOfVertices(void)=0; 174 182 175 virtual int Id()=0;176 virtual int Sid()=0;177 183 virtual bool IsNodeOnShelfFromFlags(IssmDouble* flags)=0; 178 184 virtual bool IsOnBed()=0; 179 185 virtual bool IsOnSurface()=0; 180 virtual bool IsIceInElement()=0;181 186 virtual void GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0; 182 187 virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0; -
issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
r17514 r17516 371 371 } 372 372 /*}}}*/ 373 /*FUNCTION Penta::DeepEcho{{{*/374 void Penta::DeepEcho(void){375 376 _printf_("Penta:\n");377 _printf_(" id: " << id << "\n");378 nodes[0]->DeepEcho();379 nodes[1]->DeepEcho();380 nodes[2]->DeepEcho();381 nodes[3]->DeepEcho();382 nodes[4]->DeepEcho();383 nodes[5]->DeepEcho();384 material->DeepEcho();385 matpar->DeepEcho();386 _printf_(" neighbor ids: " << verticalneighbors[0]->Id() << "-" << verticalneighbors[1]->Id() << "\n");387 _printf_(" parameters\n");388 parameters->DeepEcho();389 _printf_(" inputs\n");390 inputs->DeepEcho();391 }392 /*}}}*/393 373 /*FUNCTION Penta::Delta18oParameterization{{{*/ 394 374 void Penta::Delta18oParameterization(void){ … … 463 443 /*clean-up*/ 464 444 delete gauss; 465 }466 /*}}}*/467 /*FUNCTION Penta::Echo{{{*/468 469 void Penta::Echo(void){470 this->DeepEcho();471 }472 /*}}}*/473 /*FUNCTION Penta::ThermalToEnthalpy{{{*/474 void Penta::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){475 matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);476 }477 /*}}}*/478 /*FUNCTION Penta::EnthalpyToThermal{{{*/479 void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){480 matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);481 }482 /*}}}*/483 /*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/484 IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){485 return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);486 }487 /*}}}*/488 /*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/489 IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){490 return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);491 445 } 492 446 /*}}}*/ … … 848 802 } 849 803 /*}}}*/ 850 /*FUNCTION Penta::GetNumberOfNodesPressure;{{{*/851 int Penta::GetNumberOfNodesPressure(void){852 return this->NumberofNodesPressure();853 }854 /*}}}*/855 /*FUNCTION Penta::GetNumberOfNodesVelocity;{{{*/856 int Penta::GetNumberOfNodesVelocity(void){857 return this->NumberofNodesVelocity();858 }859 /*}}}*/860 804 /*FUNCTION Penta::GetNumberOfVertices;{{{*/ 861 805 int Penta::GetNumberOfVertices(void){ … … 1160 1104 /*Assign output pointer*/ 1161 1105 *pxyz_zero= xyz_zero; 1162 }1163 /*}}}*/1164 /*FUNCTION Penta::Sid {{{*/1165 int Penta::Sid(){1166 1167 return sid;1168 1169 }1170 /*}}}*/1171 /*FUNCTION Penta::Id {{{*/1172 int Penta::Id(void){1173 return id;1174 1106 } 1175 1107 /*}}}*/ … … 1643 1575 } 1644 1576 /*}}}*/ 1645 /*FUNCTION Penta::IsIceInElement {{{*/1646 bool Penta::IsIceInElement(){1647 1648 /*Get levelset*/1649 IssmDouble ls[NUMVERTICES];1650 GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);1651 1652 /*If the level set one one of the nodes is <0, ice is present in this element*/1653 if(ls[0]<0. || ls[1]<0. || ls[2]<0.)1654 return true;1655 else1656 return false;1657 }1658 /*}}}*/1659 1577 /*FUNCTION Penta::MinEdgeLength{{{*/ 1660 1578 IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){ … … 1813 1731 void Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ 1814 1732 1815 int i;1816 1733 IssmDouble v13[3],v23[3]; 1817 1734 IssmDouble normal[3]; 1818 1735 IssmDouble normal_norm; 1819 1736 1820 for (i=0;i<3;i++){1737 for(int i=0;i<3;i++){ 1821 1738 v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; 1822 1739 v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; … … 1826 1743 normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; 1827 1744 normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; 1828 normal_norm=sqrt( pow(normal[0],2)+pow(normal[1],2)+pow(normal[2],2));1745 normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]); 1829 1746 1830 1747 /*Bed normal is opposite to surface normal*/ … … 1917 1834 /*clean-up*/ 1918 1835 delete gauss; 1919 }1920 /*}}}*/1921 /*FUNCTION Penta::PureIceEnthalpy{{{*/1922 IssmDouble Penta::PureIceEnthalpy(IssmDouble pressure){1923 1924 return this->matpar->PureIceEnthalpy(pressure);1925 1836 } 1926 1837 /*}}}*/ … … 2070 1981 void Penta::SetTemporaryElementType(int element_type_in){ 2071 1982 this->element_type=element_type_in; 2072 }2073 /*}}}*/2074 /*FUNCTION Penta::SetwiseNodeConnectivity{{{*/2075 void Penta::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){2076 2077 /*Intermediaries*/2078 const int numnodes = this->NumberofNodes();2079 2080 /*Output */2081 int d_nz = 0;2082 int o_nz = 0;2083 2084 /*Loop over all nodes*/2085 for(int i=0;i<numnodes;i++){2086 2087 if(!flags[this->nodes[i]->Lid()]){2088 2089 /*flag current node so that no other element processes it*/2090 flags[this->nodes[i]->Lid()]=true;2091 2092 int counter=0;2093 while(flagsindices[counter]>=0) counter++;2094 flagsindices[counter]=this->nodes[i]->Lid();2095 2096 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/2097 switch(set2_enum){2098 case FsetEnum:2099 if(nodes[i]->indexing.fsize){2100 if(this->nodes[i]->IsClone())2101 o_nz += 1;2102 else2103 d_nz += 1;2104 }2105 break;2106 case GsetEnum:2107 if(nodes[i]->indexing.gsize){2108 if(this->nodes[i]->IsClone())2109 o_nz += 1;2110 else2111 d_nz += 1;2112 }2113 break;2114 case SsetEnum:2115 if(nodes[i]->indexing.ssize){2116 if(this->nodes[i]->IsClone())2117 o_nz += 1;2118 else2119 d_nz += 1;2120 }2121 break;2122 default: _error_("not supported");2123 }2124 }2125 }2126 2127 /*Special case: 2d/3d coupling, the node of this element might be connected2128 *to the basal element*/2129 int analysis_type,approximation,numlayers;2130 parameters->FindParam(&analysis_type,AnalysisTypeEnum);2131 if(analysis_type==StressbalanceAnalysisEnum){2132 inputs->GetInputValue(&approximation,ApproximationEnum);2133 if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){2134 parameters->FindParam(&numlayers,MeshNumberoflayersEnum);2135 o_nz += numlayers*3;2136 d_nz += numlayers*3;2137 }2138 }2139 2140 /*Assign output pointers: */2141 *pd_nz=d_nz;2142 *po_nz=o_nz;2143 1983 } 2144 1984 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/Penta.h
r17514 r17516 32 32 public: 33 33 34 int id;35 int sid;36 37 34 Penta **verticalneighbors; // 2 neighbors: first one under, second one above 38 35 … … 44 41 /*Object virtual functions definitions: {{{*/ 45 42 Object *copy(); 46 void DeepEcho();47 void Echo();48 43 int ObjectEnum(); 49 int Id();50 44 /*}}}*/ 51 45 /*Update virtual functions definitions: {{{*/ … … 67 61 int FiniteElement(void); 68 62 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); 69 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);70 63 void Delta18oParameterization(void); 71 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure);72 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure);73 64 Penta* GetUpperPenta(void); 74 65 Penta* GetLowerPenta(void); … … 83 74 int GetNodeIndex(Node* node); 84 75 int GetNumberOfNodes(void); 85 int GetNumberOfNodesPressure(void);86 int GetNumberOfNodesVelocity(void);87 76 int GetNumberOfVertices(void); 88 77 void GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type); … … 94 83 void GetVerticesCoordinatesTop(IssmDouble** pxyz_list); 95 84 96 int Sid();97 85 void InputDepthAverageAtBase(int enum_type,int average_enum_type); 98 86 void InputDuplicate(int original_enum,int new_enum); … … 101 89 int NumberofNodesPressure(void); 102 90 int VelocityInterpolation(); 103 IssmDouble PureIceEnthalpy(IssmDouble pressure);104 91 int PressureInterpolation(); 105 92 bool IsZeroLevelset(int levelset_enum); … … 189 176 void NormalTop(IssmDouble* bed_normal, IssmDouble* xyz_list); 190 177 ElementMatrix* CreateBasalMassMatrix(void); 191 IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure);192 IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure);193 178 void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); 194 179 … … 208 193 void JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); 209 194 void JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss); 210 bool IsIceInElement(void);211 195 Gauss* NewGauss(void); 212 196 Gauss* NewGauss(int order); -
issm/trunk-jpl/src/c/classes/Elements/Seg.cpp
r17472 r17516 65 65 } 66 66 /*}}}*/ 67 /*FUNCTION Seg::Echo{{{*/68 void Seg::Echo(void){69 _printf_("Seg:\n");70 _printf_(" id: " << id << "\n");71 if(nodes){72 nodes[0]->Echo();73 nodes[1]->Echo();74 }75 else _printf_("nodes = NULL\n");76 77 if (material) material->Echo();78 else _printf_("material = NULL\n");79 80 if (matpar) matpar->Echo();81 else _printf_("matpar = NULL\n");82 83 _printf_(" parameters\n");84 if (parameters) parameters->Echo();85 else _printf_("parameters = NULL\n");86 87 _printf_(" inputs\n");88 if (inputs) inputs->Echo();89 else _printf_("inputs=NULL\n");90 }91 /*}}}*/92 67 /*FUNCTION Seg::FiniteElement{{{*/ 93 68 int Seg::FiniteElement(void){ … … 95 70 } 96 71 /*}}}*/ 97 /*FUNCTION Seg::DeepEcho{{{*/98 void Seg::DeepEcho(void){99 100 _printf_("Seg:\n");101 _printf_(" id: " << id << "\n");102 if(nodes){103 nodes[0]->DeepEcho();104 nodes[1]->DeepEcho();105 }106 else _printf_("nodes = NULL\n");107 108 if (material) material->DeepEcho();109 else _printf_("material = NULL\n");110 111 if (matpar) matpar->DeepEcho();112 else _printf_("matpar = NULL\n");113 114 _printf_(" parameters\n");115 if (parameters) parameters->DeepEcho();116 else _printf_("parameters = NULL\n");117 118 _printf_(" inputs\n");119 if (inputs) inputs->DeepEcho();120 else _printf_("inputs=NULL\n");121 122 return;123 }124 /*}}}*/125 72 /*FUNCTION Seg::ObjectEnum{{{*/ 126 73 int Seg::ObjectEnum(void){ 127 74 128 75 return SegEnum; 129 130 }131 /*}}}*/132 /*FUNCTION Seg::Id {{{*/133 int Seg::Id(){134 135 return id;136 76 137 77 } … … 186 126 187 127 }/*}}}*/ 188 /*FUNCTION Seg::IsIceInElement {{{*/189 bool Seg::IsIceInElement(){190 191 /*Get levelset*/192 IssmDouble ls[NUMVERTICES];193 GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);194 195 /*If the level set on one of the nodes is <0, ice is present in this element*/196 if(ls[0]<0. || ls[1]<0.)197 return true;198 else199 return false;200 }201 /*}}}*/202 128 bool Seg::IsIcefront(void){/*{{{*/ 203 129 -
issm/trunk-jpl/src/c/classes/Elements/Seg.h
r17472 r17516 30 30 public: 31 31 32 int id;33 int sid;34 35 32 /*Seg constructors, destructors {{{*/ 36 33 Seg(){}; … … 39 36 /*}}}*/ 40 37 /*Object virtual functions definitions:{{{ */ 41 void Echo();42 void DeepEcho();43 int Id();44 38 int ObjectEnum(); 45 39 Object *copy(); … … 64 58 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; 65 59 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; 66 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");};67 60 void Delta18oParameterization(void){_error_("not implemented yet");}; 68 61 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; 69 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};70 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};71 IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};72 IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};73 62 int FiniteElement(void); 74 63 Element* GetUpperElement(void){_error_("not implemented yet");}; … … 78 67 int GetNodeIndex(Node* node){_error_("not implemented yet");}; 79 68 int GetNumberOfNodes(void); 80 int GetNumberOfNodesVelocity(void){_error_("not implemented yet");};81 int GetNumberOfNodesPressure(void){_error_("not implemented yet");};82 69 int GetNumberOfVertices(void); 83 70 void GetVerticesCoordinates(IssmDouble** pxyz_list); 84 71 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list){_error_("not implemented yet");}; 85 72 void GetVerticesCoordinatesTop(IssmDouble** pxyz_list){_error_("not implemented yet");}; 86 int Sid(){_error_("not implemented yet");};87 73 bool IsOnBed(){_error_("not implemented yet");}; 88 74 bool IsOnSurface(){_error_("not implemented yet");}; … … 102 88 void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; 103 89 void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; 104 bool IsIceInElement();105 90 void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); 106 91 void NormalTop(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");}; … … 111 96 Element* SpawnTopElement(void){_error_("not implemented yet");}; 112 97 IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; 113 IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};114 98 int PressureInterpolation(void){_error_("not implemented yet");}; 115 99 void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");}; -
issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
r17515 r17516 52 52 /*}}}*/ 53 53 54 /*FUNCTION Tetra::Echo{{{*/55 void Tetra::Echo(void){56 _printf_("Tetra:\n");57 _printf_(" id: " << id << "\n");58 if(nodes){59 nodes[0]->Echo();60 nodes[1]->Echo();61 nodes[2]->Echo();62 nodes[3]->Echo();63 }64 else _printf_("nodes = NULL\n");65 66 if (material) material->Echo();67 else _printf_("material = NULL\n");68 69 if (matpar) matpar->Echo();70 else _printf_("matpar = NULL\n");71 72 _printf_(" parameters\n");73 if (parameters) parameters->Echo();74 else _printf_("parameters = NULL\n");75 76 _printf_(" inputs\n");77 if (inputs) inputs->Echo();78 else _printf_("inputs=NULL\n");79 }80 /*}}}*/81 54 /*FUNCTION Tetra::FiniteElement{{{*/ 82 55 int Tetra::FiniteElement(void){ 83 56 return this->element_type; 84 } 85 /*}}}*/ 86 /*FUNCTION Tetra::DeepEcho{{{*/ 87 void Tetra::DeepEcho(void){ 88 89 _printf_("Tetra:\n"); 90 _printf_(" id: " << id << "\n"); 91 if(nodes){ 92 nodes[0]->DeepEcho(); 93 nodes[1]->DeepEcho(); 94 nodes[2]->DeepEcho(); 95 nodes[3]->DeepEcho(); 96 } 97 else _printf_("nodes = NULL\n"); 98 99 if (material) material->DeepEcho(); 100 else _printf_("material = NULL\n"); 101 102 if (matpar) matpar->DeepEcho(); 103 else _printf_("matpar = NULL\n"); 104 105 _printf_(" parameters\n"); 106 if (parameters) parameters->DeepEcho(); 107 else _printf_("parameters = NULL\n"); 108 109 _printf_(" inputs\n"); 110 if (inputs) inputs->DeepEcho(); 111 else _printf_("inputs=NULL\n"); 112 113 return; 114 } 115 /*}}}*/ 57 } /*}}}*/ 116 58 /*FUNCTION Tetra::ObjectEnum{{{*/ 117 59 int Tetra::ObjectEnum(void){ … … 119 61 return TetraEnum; 120 62 121 } 122 /*}}}*/ 123 /*FUNCTION Tetra::Id {{{*/ 124 int Tetra::Id(){ 125 return id; 126 } 127 /*}}}*/ 63 }/*}}}*/ 128 64 129 65 /*FUNCTION Tetra::AddInput{{{*/ … … 240 176 } 241 177 /*}}}*/ 242 /*FUNCTION Tetra::GetNumberOfNodesPressure THIS ONE (and corresponding TetraRef function){{{*/243 int Tetra::GetNumberOfNodesPressure(void){244 return this->NumberofNodesPressure();245 }246 /*}}}*/247 /*FUNCTION Tetra::GetNumberOfNodesVelocity;{{{*/248 int Tetra::GetNumberOfNodesVelocity(void){249 return this->NumberofNodesVelocity();250 }251 /*}}}*/252 178 /*FUNCTION Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){{{*/ 253 179 void Tetra::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){ … … 339 265 return false; 340 266 } 341 }342 /*}}}*/343 /*FUNCTION Tetra::HasNodeOnBed THIS ONE{{{*/344 bool Tetra::HasNodeOnBed(){345 346 IssmDouble values[NUMVERTICES];347 IssmDouble sum;348 349 /*Retrieve all inputs and parameters*/350 GetInputListOnVertices(&values[0],MeshVertexonbedEnum);351 sum = values[0]+values[1]+values[2]+values[3];352 353 return sum>0.;354 267 } 355 268 /*}}}*/ … … 485 398 } 486 399 /*}}}*/ 487 /*FUNCTION Tetra::IsIceInElement THIS ONE{{{*/488 bool Tetra::IsIceInElement(){489 490 /*Get levelset*/491 IssmDouble ls[NUMVERTICES];492 GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);493 494 /*If the level set on one of the nodes is <0, ice is present in this element*/495 for(int i=0;i<NUMVERTICES;i++){496 if(ls[i]<0.) return true;497 }498 499 return false;500 }501 /*}}}*/502 400 /*FUNCTION Tetra::IsOnBed {{{*/ 503 401 bool Tetra::IsOnBed(){ … … 651 549 } 652 550 /*}}}*/ 551 /*FUNCTION Tetra::NormalBase (THIS ONE){{{*/ 552 void Tetra::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){ 553 554 IssmDouble v13[3],v23[3]; 555 IssmDouble normal[3]; 556 IssmDouble normal_norm; 557 558 for(int i=0;i<3;i++){ 559 v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; 560 v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; 561 } 562 563 normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; 564 normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; 565 normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; 566 normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]); 567 568 /*Bed normal is opposite to surface normal*/ 569 bed_normal[0]=-normal[0]/normal_norm; 570 bed_normal[1]=-normal[1]/normal_norm; 571 bed_normal[2]=-normal[2]/normal_norm; 572 } 573 /*}}}*/ 574 /*FUNCTION Tetra::NormalTop (THIS ONE){{{*/ 575 void Tetra::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){ 576 577 IssmDouble v13[3],v23[3]; 578 IssmDouble normal[3]; 579 IssmDouble normal_norm; 580 581 for(int i=0;i<3;i++){ 582 v13[i]=xyz_list[0*3+i]-xyz_list[2*3+i]; 583 v23[i]=xyz_list[1*3+i]-xyz_list[2*3+i]; 584 } 585 586 normal[0]=v13[1]*v23[2]-v13[2]*v23[1]; 587 normal[1]=v13[2]*v23[0]-v13[0]*v23[2]; 588 normal[2]=v13[0]*v23[1]-v13[1]*v23[0]; 589 normal_norm=sqrt(normal[0]*normal[0]+ normal[1]*normal[1]+ normal[2]*normal[2]); 590 591 top_normal[0]=normal[0]/normal_norm; 592 top_normal[1]=normal[1]/normal_norm; 593 top_normal[2]=normal[2]/normal_norm; 594 } 595 /*}}}*/ 653 596 /*FUNCTION Tetra::NumberofNodesPressure{{{*/ 654 597 int Tetra::NumberofNodesPressure(void){ … … 666 609 if(pe){ 667 610 if(this->element_type==MINIcondensedEnum){ 668 _error_("Not implemented"); 611 int indices[3]={12,13,14}; 612 pe->StaticCondensation(Ke,3,&indices[0]); 669 613 } 670 614 else if(this->element_type==P1bubblecondensedEnum){ … … 681 625 if(Ke){ 682 626 if(this->element_type==MINIcondensedEnum){ 683 _error_("Not implemented"); 627 int indices[3]={12,13,14}; 628 Ke->StaticCondensation(3,&indices[0]); 684 629 } 685 630 else if(this->element_type==P1bubblecondensedEnum){ … … 768 713 if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp(); 769 714 else this->nodes=NULL; 770 771 }772 /*}}}*/773 /*FUNCTION Tetra::SetwiseNodeConnectivity THIS ONE (take from Penta.cpp){{{*/774 void Tetra::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){775 776 /*Intermediaries*/777 const int numnodes = this->NumberofNodes();778 779 /*Output */780 int d_nz = 0;781 int o_nz = 0;782 783 /*Loop over all nodes*/784 for(int i=0;i<numnodes;i++){785 786 if(!flags[this->nodes[i]->Lid()]){787 788 /*flag current node so that no other element processes it*/789 flags[this->nodes[i]->Lid()]=true;790 791 int counter=0;792 while(flagsindices[counter]>=0) counter++;793 flagsindices[counter]=this->nodes[i]->Lid();794 795 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/796 switch(set2_enum){797 case FsetEnum:798 if(nodes[i]->indexing.fsize){799 if(this->nodes[i]->IsClone())800 o_nz += 1;801 else802 d_nz += 1;803 }804 break;805 case GsetEnum:806 if(nodes[i]->indexing.gsize){807 if(this->nodes[i]->IsClone())808 o_nz += 1;809 else810 d_nz += 1;811 }812 break;813 case SsetEnum:814 if(nodes[i]->indexing.ssize){815 if(this->nodes[i]->IsClone())816 o_nz += 1;817 else818 d_nz += 1;819 }820 break;821 default: _error_("not supported");822 }823 }824 }825 826 /*Assign output pointers: */827 *pd_nz=d_nz;828 *po_nz=o_nz;829 }830 /*}}}*/831 /*FUNCTION Tetra::Sid THIS ONE{{{*/832 int Tetra::Sid(){833 834 return sid;835 715 836 716 } -
issm/trunk-jpl/src/c/classes/Elements/Tetra.h
r17515 r17516 30 30 public: 31 31 32 int id;33 int sid;34 35 32 /*Tetra constructors, destructors {{{*/ 36 33 Tetra(){}; … … 39 36 /*}}}*/ 40 37 /*Object virtual functions definitions:{{{ */ 41 void Echo();42 void DeepEcho();43 int Id();44 38 int ObjectEnum(); 45 39 Object *copy(); … … 64 58 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); 65 59 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); 66 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);67 60 void Delta18oParameterization(void){_error_("not implemented yet");}; 68 61 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; 69 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};70 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};71 IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};72 IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};73 62 void FaceOnFrontIndices(int* pindex1,int* pindex2,int* pindex3); 74 63 void FaceOnBedIndices(int* pindex1,int* pindex2,int* pindex3); … … 81 70 int GetNodeIndex(Node* node){_error_("not implemented yet");}; 82 71 int GetNumberOfNodes(void); 83 int GetNumberOfNodesVelocity(void);84 int GetNumberOfNodesPressure(void);85 72 int GetNumberOfVertices(void); 86 73 void GetVerticesCoordinatesBase(IssmDouble** pxyz_list); … … 88 75 bool HasFaceOnBed(); 89 76 bool HasFaceOnSurface(); 90 bool HasNodeOnBed();91 int Sid();92 77 bool IsOnBed(); 93 78 bool IsOnSurface(); … … 107 92 void NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");}; 108 93 void NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss); 109 bool IsIceInElement();110 94 void NormalSection(IssmDouble* normal,IssmDouble* xyz_list); 111 void NormalTop(IssmDouble* normal,IssmDouble* xyz_list) {_error_("not implemented yet");};112 void NormalBase(IssmDouble* normal,IssmDouble* xyz_list) {_error_("not implemented yet");};95 void NormalTop(IssmDouble* normal,IssmDouble* xyz_list); 96 void NormalBase(IssmDouble* normal,IssmDouble* xyz_list); 113 97 int NumberofNodesVelocity(void); 114 98 int NumberofNodesPressure(void); … … 117 101 Tria* SpawnTria(int index1,int index2,int index3); 118 102 IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");}; 119 IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};120 103 int PressureInterpolation(void); 121 104 void ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){_error_("not implemented yet");}; … … 228 211 /*}}}*/ 229 212 }; 230 #endif /* _ SEG_H*/213 #endif /* _TETRA_H_*/ -
issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
r17472 r17516 95 95 96 96 /*Other*/ 97 /*FUNCTION Tria::SetwiseNodeConnectivity{{{*/98 void Tria::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){99 100 /*Intermediaries*/101 const int numnodes = this->NumberofNodes();102 103 /*Output */104 int d_nz = 0;105 int o_nz = 0;106 107 /*Loop over all nodes*/108 for(int i=0;i<numnodes;i++){109 110 if(!flags[this->nodes[i]->Lid()]){111 112 /*flag current node so that no other element processes it*/113 flags[this->nodes[i]->Lid()]=true;114 115 int counter=0;116 while(flagsindices[counter]>=0) counter++;117 flagsindices[counter]=this->nodes[i]->Lid();118 119 /*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/120 switch(set2_enum){121 case FsetEnum:122 if(nodes[i]->indexing.fsize){123 if(this->nodes[i]->IsClone())124 o_nz += 1;125 else126 d_nz += 1;127 }128 break;129 case GsetEnum:130 if(nodes[i]->indexing.gsize){131 if(this->nodes[i]->IsClone())132 o_nz += 1;133 else134 d_nz += 1;135 }136 break;137 case SsetEnum:138 if(nodes[i]->indexing.ssize){139 if(this->nodes[i]->IsClone())140 o_nz += 1;141 else142 d_nz += 1;143 }144 break;145 default: _error_("not supported");146 }147 }148 }149 150 /*Assign output pointers: */151 *pd_nz=d_nz;152 *po_nz=o_nz;153 }154 /*}}}*/155 97 /*FUNCTION Tria::AddBasalInput{{{*/ 156 98 void Tria::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){ … … 292 234 } 293 235 /*}}}*/ 294 /*FUNCTION Tria::DeepEcho{{{*/295 void Tria::DeepEcho(void){296 297 _printf_("Tria:\n");298 _printf_(" id: " << id << "\n");299 if(nodes){300 nodes[0]->DeepEcho();301 nodes[1]->DeepEcho();302 nodes[2]->DeepEcho();303 }304 else _printf_("nodes = NULL\n");305 306 if (material) material->DeepEcho();307 else _printf_("material = NULL\n");308 309 if (matpar) matpar->DeepEcho();310 else _printf_("matpar = NULL\n");311 312 _printf_(" parameters\n");313 if (parameters) parameters->DeepEcho();314 else _printf_("parameters = NULL\n");315 316 _printf_(" inputs\n");317 if (inputs) inputs->DeepEcho();318 else _printf_("inputs=NULL\n");319 320 return;321 }322 /*}}}*/323 236 /*FUNCTION Tria::Delta18oParameterization{{{*/ 324 237 void Tria::Delta18oParameterization(void){ … … 414 327 *hy=ymax-ymin; 415 328 *hz=0.; 416 }417 /*}}}*/418 /*FUNCTION Tria::Echo{{{*/419 void Tria::Echo(void){420 _printf_("Tria:\n");421 _printf_(" id: " << id << "\n");422 if(nodes){423 nodes[0]->Echo();424 nodes[1]->Echo();425 nodes[2]->Echo();426 }427 else _printf_("nodes = NULL\n");428 429 if (material) material->Echo();430 else _printf_("material = NULL\n");431 432 if (matpar) matpar->Echo();433 else _printf_("matpar = NULL\n");434 435 _printf_(" parameters\n");436 if (parameters) parameters->Echo();437 else _printf_("parameters = NULL\n");438 439 _printf_(" inputs\n");440 if (inputs) inputs->Echo();441 else _printf_("inputs=NULL\n");442 329 } 443 330 /*}}}*/ … … 922 809 } 923 810 /*}}}*/ 924 /*FUNCTION Tria::GetNumberOfNodesPressure;{{{*/925 int Tria::GetNumberOfNodesPressure(void){926 return this->NumberofNodesPressure();927 }928 /*}}}*/929 /*FUNCTION Tria::GetNumberOfNodesVelocity;{{{*/930 int Tria::GetNumberOfNodesVelocity(void){931 return this->NumberofNodesVelocity();932 }933 /*}}}*/934 811 /*FUNCTION Tria::GetNumberOfVertices;{{{*/ 935 812 int Tria::GetNumberOfVertices(void){ … … 1008 885 1009 886 return y; 1010 }1011 /*}}}*/1012 /*FUNCTION Tria::Id {{{*/1013 int Tria::Id(){1014 1015 return id;1016 1017 }1018 /*}}}*/1019 /*FUNCTION Tria::Sid {{{*/1020 int Tria::Sid(){1021 1022 return sid;1023 1024 887 } 1025 888 /*}}}*/ … … 1351 1214 } 1352 1215 /*}}}*/ 1353 /*FUNCTION Tria::HasNodeOnBed {{{*/1354 bool Tria::HasNodeOnBed(){1355 1356 IssmDouble values[NUMVERTICES];1357 IssmDouble sum;1358 1359 /*Retrieve all inputs and parameters*/1360 GetInputListOnVertices(&values[0],MeshVertexonbedEnum);1361 sum = values[0]+values[1]+values[2];1362 1363 return sum>0.;1364 }1365 /*}}}*/1366 1216 /*FUNCTION Tria::HasEdgeOnSurface {{{*/ 1367 1217 bool Tria::HasEdgeOnSurface(){ … … 1384 1234 return false; 1385 1235 } 1386 }1387 /*}}}*/1388 /*FUNCTION Tria::HasNodeOnSurface {{{*/1389 bool Tria::HasNodeOnSurface(){1390 1391 IssmDouble values[NUMVERTICES];1392 IssmDouble sum;1393 1394 /*Retrieve all inputs and parameters*/1395 GetInputListOnVertices(&values[0],MeshVertexonbedEnum);1396 sum = values[0]+values[1]+values[2];1397 1398 return sum>0.;1399 1236 } 1400 1237 /*}}}*/ … … 1540 1377 this->EdgeOnSurfaceIndices(&indices[0],&indices[1]); 1541 1378 return new GaussTria(indices[0],indices[1],order); 1542 }1543 /*}}}*/1544 /*FUNCTION Tria::IsIceInElement {{{*/1545 bool Tria::IsIceInElement(){1546 1547 /*Get levelset*/1548 IssmDouble ls[NUMVERTICES];1549 GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);1550 1551 /*If the level set on one of the nodes is <0, ice is present in this element*/1552 if(ls[0]<0. || ls[1]<0. || ls[2]<0.)1553 return true;1554 else1555 return false;1556 1379 } 1557 1380 /*}}}*/ -
issm/trunk-jpl/src/c/classes/Elements/Tria.h
r17472 r17516 32 32 public: 33 33 34 int id;35 int sid;36 37 34 /*Tria constructors, destructors {{{*/ 38 35 Tria(){}; … … 41 38 /*}}}*/ 42 39 /*Object virtual functions definitions:{{{ */ 43 void Echo();44 void DeepEcho();45 int Id();46 40 int ObjectEnum(); 47 41 Object *copy(); … … 63 57 void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); 64 58 void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); 65 void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum);66 59 void Delta18oParameterization(void); 67 60 void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); 68 void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");};69 void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};70 61 int FiniteElement(void); 71 62 Element* GetUpperElement(void){_error_("not implemented yet");}; … … 77 68 int GetNodeIndex(Node* node); 78 69 int GetNumberOfNodes(void); 79 int GetNumberOfNodesPressure(void);80 int GetNumberOfNodesVelocity(void);81 70 int GetNumberOfVertices(void); 82 int Sid();83 71 bool IsOnBed(); 84 72 bool IsOnSurface(); 85 73 bool HasEdgeOnBed(); 86 bool HasNodeOnBed();87 74 bool HasEdgeOnSurface(); 88 bool HasNodeOnSurface();89 75 void EdgeOnSurfaceIndices(int* pindex1,int* pindex); 90 76 void EdgeOnBedIndices(int* pindex1,int* pindex); … … 94 80 int NumberofNodesVelocity(void); 95 81 int NumberofNodesPressure(void); 96 bool IsIceInElement();97 82 void GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type); 98 83 void GetVectorFromInputs(Vector<IssmDouble>* vector, int name_enum); … … 111 96 Element* SpawnTopElement(void); 112 97 int VelocityInterpolation(); 113 IssmDouble PureIceEnthalpy(IssmDouble pressure){_error_("not implemented yet");};114 98 int PressureInterpolation(); 115 99 IssmDouble SurfaceArea(void); … … 192 176 void AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum); 193 177 void AddInput(int input_enum, IssmDouble* values, int interpolation_enum); 194 IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");};195 IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");};196 178 IssmDouble GetArea(void); 197 179 void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints); -
issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp
r17444 r17516 242 242 case HydrologyDCInefficientAnalysisEnum: 243 243 Ke=PenaltyCreateKMatrixHydrologyDCInefficient(kmax); 244 break;245 case DamageEvolutionAnalysisEnum:246 Ke=PenaltyCreateKMatrixDamageEvolution(kmax);247 244 break; 248 245 default: … … 276 273 case HydrologyDCInefficientAnalysisEnum: 277 274 pe=PenaltyCreatePVectorHydrologyDCInefficient(kmax); 278 break;279 case DamageEvolutionAnalysisEnum:280 pe=PenaltyCreatePVectorDamageEvolution(kmax);281 275 break; 282 276 default: … … 417 411 return; 418 412 } 419 else if (analysis_type==DamageEvolutionAnalysisEnum){420 ConstraintActivateDamageEvolution(punstable);421 return;422 }423 424 413 else{ 425 414 _error_("analysis: " << EnumToStringx(analysis_type) << " not supported yet"); … … 710 699 } 711 700 /*}}}*/ 712 /*FUNCTION Pengrid::ConstraintActivateDamageEvolution {{{*/713 void Pengrid::ConstraintActivateDamageEvolution(int* punstable){714 715 // The penalty is stable if it doesn't change during to successive iterations.716 IssmDouble max_damage;717 IssmDouble damage;718 int new_active;719 int unstable=0;720 int penalty_lock;721 722 /*check that pengrid is not a clone (penalty to be added only once)*/723 if (node->IsClone()){724 unstable=0;725 *punstable=unstable;726 return;727 }728 729 //First recover damage using the element: */730 element->GetInputValue(&damage,node,DamageDEnum);731 732 //Recover our data:733 parameters->FindParam(&penalty_lock,DamagePenaltyLockEnum);734 parameters->FindParam(&max_damage,DamageMaxDamageEnum);735 736 //Figure out if damage>max_damage, in which case, this penalty needs to be activated.737 //Would need to do the same for damage<0 if penalties are used. For now, ConstraintStatex738 //is not called in solutionsequence_damage_nonlinear, so no penalties are applied.739 740 if (damage>max_damage){741 new_active=1;742 }743 else{744 new_active=0;745 }746 747 //Figure out stability of this penalty748 if (active==new_active){749 unstable=0;750 }751 else{752 unstable=1;753 if(penalty_lock)zigzag_counter++;754 }755 756 /*If penalty keeps zigzagging more than penalty_lock times: */757 if(penalty_lock){758 if(zigzag_counter>penalty_lock){759 unstable=0;760 active=1;761 }762 }763 764 //Set penalty flag765 active=new_active;766 767 //*Assign output pointers:*/768 *punstable=unstable;769 }770 /*}}}*/771 /*FUNCTION Pengrid::PenaltyCreateKMatrixDamageEvolution {{{*/772 ElementMatrix* Pengrid::PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax){773 774 IssmDouble penalty_factor;775 776 /*Initialize Element matrix and return if necessary*/777 if(!this->active) return NULL;778 ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters);779 780 /*recover parameters: */781 parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);782 783 Ke->values[0]=kmax*pow(10.,penalty_factor);784 785 /*Clean up and return*/786 return Ke;787 }788 /*}}}*/789 /*FUNCTION Pengrid::PenaltyCreatePVectorDamageEvolution {{{*/790 ElementVector* Pengrid::PenaltyCreatePVectorDamageEvolution(IssmDouble kmax){791 792 IssmDouble penalty_factor;793 IssmDouble max_damage;794 795 /*Initialize Element matrix and return if necessary*/796 if(!this->active) return NULL;797 ElementVector* pe=new ElementVector(&node,1,this->parameters);798 799 //Recover our data:800 parameters->FindParam(&penalty_factor,DamagePenaltyFactorEnum);801 parameters->FindParam(&max_damage,DamageMaxDamageEnum);802 803 //right hand side penalizes to max_damage804 pe->values[0]=kmax*pow(10.,penalty_factor)*max_damage;805 806 /*Clean up and return*/807 return pe;808 }809 /*}}}*/810 701 /*FUNCTION Pengrid::CreatePVectorHydrologyDCInefficient {{{*/ 811 702 ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){ -
issm/trunk-jpl/src/c/classes/Loads/Pengrid.h
r17444 r17516 87 87 ElementVector* PenaltyCreatePVectorMelting(IssmDouble kmax); 88 88 void ConstraintActivateThermal(int* punstable); 89 ElementMatrix* PenaltyCreateKMatrixDamageEvolution(IssmDouble kmax);90 ElementVector* PenaltyCreatePVectorDamageEvolution(IssmDouble kmax);91 void ConstraintActivateDamageEvolution(int* punstable);92 89 ElementMatrix* PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax); 93 90 ElementVector* PenaltyCreatePVectorHydrologyDCInefficient(IssmDouble kmax); -
issm/trunk-jpl/src/c/solutionsequences/solutionsequences.h
r16188 r17516 13 13 14 14 void solutionsequence_thermal_nonlinear(FemModel* femmodel); 15 void solutionsequence_damage_nonlinear(FemModel* femmodel);16 15 void solutionsequence_hydro_nonlinear(FemModel* femmodel); 17 16 void solutionsequence_nonlinear(FemModel* femmodel,bool conserve_loads);
Note:
See TracChangeset
for help on using the changeset viewer.
