Changeset 7089
- Timestamp:
- 01/14/11 08:14:31 (14 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 18 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/Container/Elements.cpp
r6850 r7089 204 204 } 205 205 /*}}}*/ 206 /*FUNCTION Elements::NumberOfElements{{{1*/ 207 int Elements::NumberOfElements(void){ 208 209 int local_nelem=0; 210 int numberofelements; 211 212 #ifdef _PARALLEL_ 213 local_nelem=this->Size(); 214 MPI_Allreduce ( (void*)&local_nelem,(void*)&numberofelements,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); 215 #else 216 numberofelements=this->Size(); 217 #endif 218 219 return numberofelements; 220 } 221 /*}}}*/ -
issm/trunk/src/c/Container/Elements.h
r6411 r7089 33 33 void ToResults(Results* results,Parameters* parameters,int step, double time); 34 34 Patch* ResultsToPatch(void); 35 int NumberOfElements(void); 35 36 /*}}}*/ 36 37 -
issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp
r6967 r7089 36 36 /*Fetch data needed: */ 37 37 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 38 IoModelFetchData(&iomodel->elementconnectivity,NULL,NULL,iomodel_handle,"elementconnectivity"); 38 39 IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements"); 39 40 IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements"); … … 48 49 49 50 /*Create and add tria element to elements dataset: */ 50 if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,i omodel,nummodels));51 else elements->AddObject(new Penta(i+1,i,i omodel,nummodels));51 if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,i,iomodel,nummodels)); 52 else elements->AddObject(new Penta(i+1,i,i,iomodel,nummodels)); 52 53 53 54 /*Create and add material property to materials dataset: */ … … 58 59 /*Free data: */ 59 60 xfree((void**)&iomodel->elements); 61 xfree((void**)&iomodel->elementconnectivity); 60 62 xfree((void**)&iomodel->upperelements); 61 63 xfree((void**)&iomodel->lowerelements); -
issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
r6946 r7089 75 75 parameters->AddObject(new BoolParam(KffEnum,iomodel->kff)); 76 76 parameters->AddObject(new BoolParam(IoGatherEnum,iomodel->io_gather)); 77 parameters->AddObject(new BoolParam(GroundingLineMigrationEnum,iomodel->grounding_line_migration)); 77 78 78 79 /*Deal with more complex parameters*/ -
issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp
r6972 r7089 39 39 IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); 40 40 IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); 41 if(iomodel->grounding_line_migration) IoModelFetchData(&iomodel->bathymetry,NULL,NULL,iomodel_handle,"bathymetry"); 41 42 if (iomodel->dim==3){ 42 43 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); -
issm/trunk/src/c/objects/Elements/Element.h
r6953 r7089 32 32 virtual void GetSolutionFromInputs(Vec solution)=0; 33 33 virtual int GetNodeIndex(Node* node)=0; 34 virtual int Sid()=0; 34 35 virtual bool IsOnShelf()=0; 36 virtual bool IsNodeOnShelf()=0; 35 37 virtual bool IsOnBed()=0; 36 38 virtual void GetParameterListOnVertices(double* pvalue,int enumtype)=0; … … 79 81 virtual bool InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums)=0; 80 82 virtual void AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part)=0; 83 virtual int* GetHorizontalNeighboorSids(void)=0; 81 84 virtual double TimeAdapt()=0; 85 virtual void MigrateGroundingLine()=0; 86 virtual void UpdateShelfStatus()=0; 82 87 83 88 /*Implementation: */ -
issm/trunk/src/c/objects/Elements/Penta.cpp
r7073 r7089 26 26 /*FUNCTION Penta::Penta(){{{1*/ 27 27 Penta::Penta(){ 28 29 int i; 30 28 31 this->nodes=NULL; 29 32 this->matice=NULL; 30 33 this->matpar=NULL; 31 this-> neighbors=NULL;34 this->verticalneighbors=NULL; 32 35 this->inputs=NULL; 33 36 this->parameters=NULL; 34 37 this->results=NULL; 38 for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF; 35 39 } 36 40 /*}}}*/ … … 43 47 /*}}}*/ 44 48 /*FUNCTION Penta::Penta(int id, int index, IoModel* iomodel,int nummodels) {{{1*/ 45 Penta::Penta(int penta_id, int index, IoModel* iomodel,int nummodels)49 Penta::Penta(int penta_id, int penta_sid, int index, IoModel* iomodel,int nummodels) 46 50 :PentaRef(nummodels) 47 51 ,PentaHook(nummodels,index+1,iomodel->numberofelements+1) //index+1: matice id, iomodel->numberofelements+1: matpar id … … 59 63 /*id: */ 60 64 this->id=penta_id; 65 this->sid=penta_sid; 61 66 62 67 /*Build neighbors list*/ … … 66 71 else penta_elements_ids[0]=(int)(iomodel->lowerelements[index]); 67 72 this->InitHookNeighbors(penta_elements_ids); 73 74 /*Build horizontalneighborsids list: */ 75 for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1; 68 76 69 77 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. … … 78 86 this->matice=NULL; 79 87 this->matpar=NULL; 80 this-> neighbors=NULL;88 this->verticalneighbors=NULL; 81 89 } 82 90 /*}}}*/ … … 104 112 /*deal with Penta copy fields: */ 105 113 penta->id=this->id; 114 penta->sid=this->sid; 106 115 if(this->inputs){ 107 116 penta->inputs=(Inputs*)this->inputs->Copy(); … … 124 133 penta->matice=(Matice*)penta->hmatice->delivers(); 125 134 penta->matpar=(Matpar*)penta->hmatpar->delivers(); 126 penta->neighbors=(Penta**)penta->hneighbors->deliverp(); 135 penta->verticalneighbors=(Penta**)penta->hneighbors->deliverp(); 136 137 /*neighbors: */ 138 for(i=0;i<3;i++)penta->horizontalneighborsids[i]=this->horizontalneighborsids[i]; 127 139 128 140 return penta; … … 154 166 /*marshall Penta data: */ 155 167 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); 168 memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid); 156 169 memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); 157 170 … … 195 208 xfree((void**)&marshalled_results); 196 209 210 /*marshall horizontal neighbors: */ 211 memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int); 212 197 213 *pmarshalled_dataset=marshalled_dataset; 198 214 return; … … 211 227 212 228 return sizeof(id) 229 +sizeof(sid) 213 230 +hnodes_size 214 231 +sizeof(numanalyses) … … 219 236 +inputs->MarshallSize() 220 237 +results->MarshallSize() 238 +3*sizeof(int) 221 239 +sizeof(int); //sizeof(int) for enum type 222 240 } … … 235 253 *object data (thanks to DataSet::Demarshall):*/ 236 254 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); 255 memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid); 237 256 memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); 238 257 … … 260 279 matice=NULL; 261 280 matpar=NULL; 262 neighbors=NULL;281 verticalneighbors=NULL; 263 282 264 283 /*demarshall inputs and results: */ … … 268 287 /*parameters: may not exist even yet, so let Configure handle it: */ 269 288 this->parameters=NULL; 289 290 /*neighbors: */ 291 memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int); 270 292 271 293 /*return: */ … … 425 447 this->matice=(Matice*)this->hmatice->delivers(); 426 448 this->matpar=(Matpar*)this->hmatpar->delivers(); 427 this-> neighbors=(Penta**)this->hneighbors->deliverp();449 this->verticalneighbors=(Penta**)this->hneighbors->deliverp(); 428 450 429 451 /*point parameters to real dataset: */ … … 499 521 500 522 /*Checks in debugging {{{2*/ 501 _assert_(this->nodes && this->matice && this->matpar && this-> neighbors && this->parameters && this->inputs);523 _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs); 502 524 /*}}}*/ 503 525 … … 1794 1816 1795 1817 /*if debugging mode, check that all pointers exist {{{2*/ 1796 _assert_(this->nodes && this->matice && this->matpar && this-> neighbors && this->parameters && this->inputs);1818 _assert_(this->nodes && this->matice && this->matpar && this->verticalneighbors && this->parameters && this->inputs); 1797 1819 /*}}}*/ 1798 1820 … … 2878 2900 matice->DeepEcho(); 2879 2901 matpar->DeepEcho(); 2880 printf(" neighbor ids: %i-%i\n", neighbors[0]->Id(),neighbors[1]->Id());2902 printf(" neighbor ids: %i-%i\n",verticalneighbors[0]->Id(),verticalneighbors[1]->Id()); 2881 2903 printf(" parameters\n"); 2882 2904 parameters->DeepEcho(); … … 2885 2907 printf(" results\n"); 2886 2908 results->DeepEcho(); 2909 printf("neighboor sids: \n"); 2910 printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]); 2911 2887 2912 return; 2888 2913 } … … 3063 3088 } 3064 3089 /*}}}*/ 3090 /*FUNCTION Penta::GetHorizontalNeighboorSids {{{1*/ 3091 int* Penta::GetHorizontalNeighboorSids(){ 3092 3093 /*return PentaRef field*/ 3094 return &this->horizontalneighborsids[0]; 3095 3096 } 3097 /*}}}*/ 3065 3098 /*FUNCTION Penta::GetLowerElement{{{1*/ 3066 3099 Penta* Penta::GetLowerElement(void){ … … 3068 3101 Penta* upper_penta=NULL; 3069 3102 3070 upper_penta=(Penta*) neighbors[0]; //first one (0) under, second one (1) above3103 upper_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above 3071 3104 3072 3105 return upper_penta; … … 3497 3530 Penta* upper_penta=NULL; 3498 3531 3499 upper_penta=(Penta*) neighbors[1]; //first one under, second one above3532 upper_penta=(Penta*)verticalneighbors[1]; //first one under, second one above 3500 3533 3501 3534 return upper_penta; … … 3533 3566 3534 3567 return z; 3568 } 3569 /*}}}*/ 3570 /*FUNCTION Penta::Sid {{{1*/ 3571 int Penta::Sid(){ 3572 3573 return sid; 3574 3535 3575 } 3536 3576 /*}}}*/ … … 5191 5231 } 5192 5232 /*}}}*/ 5233 /*FUNCTION Penta::IsNodeOnShelf {{{1*/ 5234 bool Penta::IsNodeOnShelf(){ 5235 5236 int i; 5237 bool shelf=false; 5238 5239 for(i=0;i<6;i++){ 5240 if (nodes[i]->IsOnShelf()){ 5241 shelf=true; 5242 break; 5243 } 5244 } 5245 return shelf; 5246 } 5247 /*}}}*/ 5193 5248 /*FUNCTION Penta::IsOnSurface{{{1*/ 5194 5249 bool Penta::IsOnSurface(void){ … … 5302 5357 /*Assign output pointers:*/ 5303 5358 *pmaxvz=maxvz; 5359 } 5360 /*}}}*/ 5361 /*FUNCTION Penta::MigrateGroundingLine{{{1*/ 5362 void Penta::MigrateGroundingLine(void){ 5363 _error_("not supported yet!"); 5304 5364 } 5305 5365 /*}}}*/ … … 5417 5477 int numrows = 0; 5418 5478 int numnodes = 0; 5479 int temp_numnodes=0; 5419 5480 5420 5481 /*Go through all the results objects, and update the counters: */ … … 5424 5485 numrows++; 5425 5486 /*now, how many vertices and how many nodal values for this result? :*/ 5426 numnodes=elementresult->NumberOfNodalValues(); //ask result object. 5487 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object. 5488 if(temp_numnodes>numnodes)numnodes=temp_numnodes; 5427 5489 } 5428 5490 … … 6071 6133 } 6072 6134 /*}}}*/ 6135 /*FUNCTION Penta::UpdateShelfStatus{{{1*/ 6136 void Penta::UpdateShelfStatus(void){ 6137 _error_("Not implemented yet"); 6138 } 6139 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Penta.h
r6987 r7089 32 32 33 33 int id; 34 int sid; 34 35 35 36 Node **nodes; // 6 nodes 36 37 Matice *matice; // 1 material ice 37 38 Matpar *matpar; // 1 material parameter 38 Penta **neighbors; // 2 neighbors: first one under, second one above 39 Penta **verticalneighbors; // 2 neighbors: first one under, second one above 40 int horizontalneighborsids[3]; 39 41 40 42 Parameters *parameters; //pointer to solution parameters … … 44 46 /*Penta constructors and destructor: {{{1*/ 45 47 Penta(); 46 Penta(int penta_id,int i, IoModel* iomodel,int nummodels);48 Penta(int penta_id,int penta_sid,int i, IoModel* iomodel,int nummodels); 47 49 ~Penta(); 48 50 /*}}}*/ … … 88 90 void GradjB(Vec gradient); 89 91 void GradjDrag(Vec gradient); 92 int Sid(); 90 93 void InputControlUpdate(double scalar,bool save_parameter); 91 94 void InputArtificialNoise(int enum_type,double min, double max); … … 106 109 void MaxVy(double* pmaxvy, bool process_units); 107 110 void MaxVz(double* pmaxvz, bool process_units); 111 void MigrateGroundingLine(); 108 112 void MinVel(double* pminvel, bool process_units); 109 113 void MinVx(double* pminvx, bool process_units); … … 121 125 double SurfaceArea(void); 122 126 void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); 127 void UpdateShelfStatus(void); 123 128 double TimeAdapt(); 129 int* GetHorizontalNeighboorSids(void); 124 130 /*}}}*/ 125 131 /*Penta specific routines:{{{1*/ … … 229 235 bool IsOnBed(void); 230 236 bool IsOnShelf(void); 237 bool IsNodeOnShelf(void); 231 238 bool IsOnWater(void); 232 239 double MinEdgeLength(double xyz_list[6][3]); -
issm/trunk/src/c/objects/Elements/Tria.cpp
r6971 r7089 26 26 Tria::Tria(){ 27 27 28 int i; 29 28 30 this->nodes=NULL; 29 31 this->matice=NULL; 30 32 this->matpar=NULL; 33 for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF; 31 34 this->inputs=NULL; 32 35 this->parameters=NULL; … … 35 38 } 36 39 /*}}}*/ 37 /*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel,int nummodels){{{1*/38 Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels)40 /*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{1*/ 41 Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels) 39 42 :TriaRef(nummodels) 40 43 ,TriaHook(nummodels,index+1,iomodel->numberofelements+1){ 41 /*id: */ 42 this->id=tria_id; 43 44 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. 45 this->parameters=NULL; 46 47 /*intialize inputs and results: */ 48 this->inputs=new Inputs(); 49 this->results=new Results(); 50 51 /*initialize pointers:*/ 52 this->nodes=NULL; 53 this->matice=NULL; 54 this->matpar=NULL; 44 45 int i; 46 /*id: */ 47 this->id=tria_id; 48 this->sid=tria_sid; 49 50 //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. 51 this->parameters=NULL; 52 53 /*Build horizontalneighborsids list: */ 54 for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1; 55 56 /*intialize inputs and results: */ 57 this->inputs=new Inputs(); 58 this->results=new Results(); 59 60 /*initialize pointers:*/ 61 this->nodes=NULL; 62 this->matice=NULL; 63 this->matpar=NULL; 55 64 56 65 } … … 84 93 /*deal with Tria fields: */ 85 94 tria->id=this->id; 95 tria->sid=this->sid; 86 96 if(this->inputs){ 87 97 tria->inputs=(Inputs*)this->inputs->Copy(); … … 104 114 tria->matice=(Matice*)tria->hmatice->delivers(); 105 115 tria->matpar=(Matpar*)tria->hmatpar->delivers(); 116 117 /*neighbors: */ 118 for(i=0;i<3;i++)tria->horizontalneighborsids[i]=this->horizontalneighborsids[i]; 106 119 107 120 return tria; … … 133 146 /*marshall Tria data: */ 134 147 memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); 148 memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid); 135 149 memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); 136 150 … … 174 188 xfree((void**)&marshalled_results); 175 189 190 /*marshall horizontal neighbors: */ 191 memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int); 192 176 193 *pmarshalled_dataset=marshalled_dataset; 177 194 return; … … 190 207 191 208 return sizeof(id) 209 +sizeof(sid) 192 210 +hnodes_size 193 211 +sizeof(numanalyses) … … 197 215 +inputs->MarshallSize() 198 216 +results->MarshallSize() 217 +3*sizeof(int) 199 218 +sizeof(int); //sizeof(int) for enum type 200 219 } … … 213 232 *object data (thanks to DataSet::Demarshall):*/ 214 233 memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); 234 memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid); 215 235 memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses); 216 236 … … 244 264 /*parameters: may not exist even yet, so let Configure handle it: */ 245 265 this->parameters=NULL; 266 267 /*neighbors: */ 268 memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int); 246 269 247 270 /*return: */ … … 2614 2637 if (results) results->DeepEcho(); 2615 2638 else printf("results=NULL\n"); 2639 2640 printf("neighboor sids: \n"); 2641 printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]); 2616 2642 2617 2643 return; … … 2654 2680 if (results) results->Echo(); 2655 2681 else printf("results=NULL\n"); 2682 2683 printf("neighboor sids: \n"); 2684 printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]); 2656 2685 } 2657 2686 /*}}}*/ … … 2715 2744 /*return TriaRef field*/ 2716 2745 return this->element_type; 2746 2747 } 2748 /*}}}*/ 2749 /*FUNCTION Tria::GetHorizontalNeighboorSids {{{1*/ 2750 int* Tria::GetHorizontalNeighboorSids(){ 2751 2752 /*return TriaRef field*/ 2753 return &this->horizontalneighborsids[0]; 2717 2754 2718 2755 } … … 3340 3377 } 3341 3378 /*}}}*/ 3379 /*FUNCTION Tria::Sid {{{1*/ 3380 int Tria::Sid(){ 3381 3382 return sid; 3383 3384 } 3385 /*}}}*/ 3342 3386 /*FUNCTION Tria::InputArtificialNoise{{{1*/ 3343 3387 void Tria::InputArtificialNoise(int enum_type,double min,double max){ … … 3548 3592 for(i=0;i<3;i++)nodeinputs[i]=iomodel->bed[tria_vertex_ids[i]-1]; 3549 3593 this->inputs->AddInput(new TriaVertexInput(BedEnum,nodeinputs)); 3594 } 3595 if(iomodel->grounding_line_migration){ 3596 for(i=0;i<3;i++)nodeinputs[i]=iomodel->bathymetry[tria_vertex_ids[i]-1]; 3597 this->inputs->AddInput(new TriaVertexInput(BathymetryEnum,nodeinputs)); 3550 3598 } 3551 3599 if (iomodel->drag_coefficient) { … … 4172 4220 } 4173 4221 /*}}}*/ 4222 /*FUNCTION Tria::IsNodeOnShelf {{{1*/ 4223 bool Tria::IsNodeOnShelf(){ 4224 4225 int i; 4226 bool shelf=false; 4227 4228 for(i=0;i<3;i++){ 4229 if (nodes[i]->IsOnShelf()){ 4230 shelf=true; 4231 break; 4232 } 4233 } 4234 return shelf; 4235 } 4236 /*}}}*/ 4174 4237 /*FUNCTION Tria::IsOnWater {{{1*/ 4175 4238 bool Tria::IsOnWater(){ … … 4335 4398 } 4336 4399 /*}}}*/ 4400 /*FUNCTION Tria::MigrateGroundingLine{{{1*/ 4401 void Tria::MigrateGroundingLine(void){ 4402 4403 int i; 4404 4405 double h[3]; 4406 double s[3]; 4407 double b[3]; 4408 double ba[3]; 4409 double bed_hydro; 4410 int isonshelf[3]; 4411 Input *surface_input = NULL; 4412 Input *bed_input = NULL; 4413 double rho_water; 4414 double rho_ice; 4415 double density; 4416 double *values = NULL; 4417 bool elementonshelf = false; 4418 4419 4420 /*Recover info at the vertices: */ 4421 surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); 4422 bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); 4423 if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); 4424 4425 GetParameterListOnVertices(&h[0],ThicknessEnum); 4426 GetParameterListOnVertices(&s[0],SurfaceEnum); 4427 GetParameterListOnVertices(&b[0],BedEnum); 4428 GetParameterListOnVertices(&ba[0],BathymetryEnum); 4429 for(i=0;i<3;i++)isonshelf[i]=nodes[i]->IsOnShelf(); 4430 4431 /*material parameters: */ 4432 rho_water=matpar->GetRhoWater(); 4433 rho_ice=matpar->GetRhoIce(); 4434 density=rho_ice/rho_water; 4435 4436 /*go through vertices, and update inputs, considering them to be TriaVertex type: */ 4437 for(i=0;i<3;i++){ 4438 if (isonshelf[i]){ 4439 /*This node is on the shelf. See if its bed is going under the bathymetry: */ 4440 if(b[i]<=ba[i]){ //<= because Neff being 0 when b=ba, drag will be 0 anyway. 4441 /*The ice shelf is getting grounded, the thickness is the same, so just update the bed to stick to the bathymetry and elevate the surface accordingly: */ 4442 b[i]=ba[i]; 4443 s[i]=b[i]+h[i]; 4444 } 4445 else{ 4446 /*do nothing, we are still floating.*/ 4447 } 4448 } 4449 else{ 4450 /*This node is on the sheet, near the grounding line. See if wants to unground. To 4451 * do so, we compute the hydrostatic bed, and if it is > bathymetry, then we unground: */ 4452 bed_hydro=-density*h[i]; 4453 if (bed_hydro>ba[i]){ 4454 /*We are now floating, bed and surface are determined from hydrostatic equilibrium: */ 4455 s[i]=(1-density)*h[i]; 4456 b[i]=-density*h[i]; 4457 } 4458 else{ 4459 /*do nothing, we are still grounded.*/ 4460 } 4461 } 4462 } 4463 4464 /*Surface and bed are updated. Update inputs:*/ 4465 surface_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=s[i]; 4466 bed_input->GetValuesPtr(&values,NULL); for(i=0;i<3;i++)values[i]=b[i]; 4467 4468 } 4469 /*}}}*/ 4337 4470 /*FUNCTION Tria::MinVel{{{1*/ 4338 4471 void Tria::MinVel(double* pminvel, bool process_units){ … … 4426 4559 int numrows = 0; 4427 4560 int numnodes = 0; 4561 int temp_numnodes = 0; 4428 4562 4429 4563 /*Go through all the results objects, and update the counters: */ … … 4433 4567 numrows++; 4434 4568 /*now, how many vertices and how many nodal values for this result? :*/ 4435 numnodes=elementresult->NumberOfNodalValues(); //ask result object. 4569 temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object. 4570 if(temp_numnodes>numnodes)numnodes=temp_numnodes; 4436 4571 } 4437 4572 … … 5196 5331 } 5197 5332 /*}}}*/ 5333 /*FUNCTION Tria::UpdateShelfStatus{{{1*/ 5334 void Tria::UpdateShelfStatus(void){ 5335 5336 int i; 5337 5338 double h[3]; 5339 double s[3]; 5340 double b[3]; 5341 double ba[3]; 5342 Input *surface_input = NULL; 5343 Input *bed_input = NULL; 5344 double rho_water; 5345 double rho_ice; 5346 double density; 5347 bool elementonshelf = false; 5348 5349 5350 /*Recover info at the vertices: */ 5351 surface_input =inputs->GetInput(SurfaceEnum); _assert_(surface_input); 5352 bed_input =inputs->GetInput(BedEnum); _assert_(bed_input); 5353 if((surface_input->Enum()!=TriaVertexInputEnum) | (bed_input->Enum()!=TriaVertexInputEnum))_error_(" not supported yet for bed and surface interpolations not P1!"); 5354 5355 GetParameterListOnVertices(&h[0],ThicknessEnum); 5356 GetParameterListOnVertices(&s[0],SurfaceEnum); 5357 GetParameterListOnVertices(&b[0],BedEnum); 5358 GetParameterListOnVertices(&ba[0],BathymetryEnum); 5359 5360 /*material parameters: */ 5361 rho_water=matpar->GetRhoWater(); 5362 rho_ice=matpar->GetRhoIce(); 5363 density=rho_ice/rho_water; 5364 5365 /*go through vertices, and figure out if they are grounded or not, then update their status: */ 5366 for(i=0;i<3;i++){ 5367 if(b[i]<=ba[i]){ 5368 nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,false)); 5369 nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,true)); 5370 } 5371 else{ 5372 nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceShelfEnum,true)); 5373 nodes[i]->inputs->AddInput(new BoolInput(NodeOnIceSheetEnum,false)); 5374 } 5375 } 5376 5377 /*Now, update shelf status of element. An element can only be on shelf if all its nodes are on shelf: */ 5378 elementonshelf=true; 5379 for(i=0;i<3;i++){ 5380 if(nodes[i]->IsOnSheet()){ 5381 elementonshelf=false; 5382 break; 5383 } 5384 } 5385 this->inputs->AddInput(new BoolInput(ElementOnIceShelfEnum,elementonshelf)); 5386 } 5387 /*}}}*/ -
issm/trunk/src/c/objects/Elements/Tria.h
r6953 r7089 30 30 31 31 int id; 32 int sid; 32 33 33 34 Node **nodes; // 3 nodes 34 35 Matice *matice; // 1 material ice 35 36 Matpar *matpar; // 1 material parameter 37 int horizontalneighborsids[3]; 36 38 37 39 Parameters *parameters; //pointer to solution parameters … … 41 43 /*Tria constructors, destructors {{{1*/ 42 44 Tria(); 43 Tria(int tria_id,int i, IoModel* iomodel,int nummodels);45 Tria(int tria_id,int tria_sid,int i, IoModel* iomodel,int nummodels); 44 46 ~Tria(); 45 47 /*}}}*/ … … 78 80 void CreatePVector(Vec pg, Vec pf); 79 81 int GetNodeIndex(Node* node); 82 int Sid(); 80 83 bool IsOnBed(); 81 84 bool IsOnShelf(); 85 bool IsNodeOnShelf(); 82 86 bool IsOnWater(); 83 87 void GetSolutionFromInputs(Vec solution); … … 109 113 void MaxVy(double* pmaxvy, bool process_units); 110 114 void MaxVz(double* pmaxvz, bool process_units); 115 void MigrateGroundingLine(); 111 116 void MinVel(double* pminvel, bool process_units); 112 117 void MinVx(double* pminvx, bool process_units); … … 124 129 double SurfaceArea(void); 125 130 void Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type); 131 void UpdateShelfStatus(void); 126 132 double TimeAdapt(); 133 int* GetHorizontalNeighboorSids(void); 127 134 /*}}}*/ 128 135 /*Tria specific routines:{{{1*/ -
issm/trunk/src/c/objects/Inputs/BoolInput.cpp
r6412 r7089 148 148 /*FUNCTION BoolInput::SpawnResult{{{1*/ 149 149 ElementResult* BoolInput::SpawnResult(int step, double time){ 150 151 _error_(" not supported yet!");150 151 return new BoolElementResult(this->enum_type,this->value,step,time); 152 152 153 153 } -
issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
r6413 r7089 399 399 400 400 *pvalues=this->values; 401 *pnum_values=3;402 403 } 404 /*}}}*/ 401 if(pnum_values)*pnum_values=3; 402 403 } 404 /*}}}*/ -
issm/trunk/src/c/objects/IoModel.cpp
r6946 r7089 38 38 xfree((void**)&this->z); 39 39 xfree((void**)&this->elements); 40 xfree((void**)&this->elementconnectivity); 40 41 xfree((void**)&this->elements_type); 41 42 xfree((void**)&this->vertices_type); … … 59 60 xfree((void**)&this->surface); 60 61 xfree((void**)&this->bed); 62 xfree((void**)&this->bathymetry); 61 63 xfree((void**)&this->vx_obs); 62 64 xfree((void**)&this->vy_obs); … … 191 193 IoModelFetchData(&this->waitonlock,iomodel_handle,"waitonlock"); 192 194 IoModelFetchData(&this->kff,iomodel_handle,"kff"); 195 IoModelFetchData(&this->grounding_line_migration,iomodel_handle,"grounding_line_migration"); 193 196 194 197 /*!Get thermal parameters: */ … … 214 217 IoModelFetchData(&this->qmu_save_femmodel,iomodel_handle,"qmu_save_femmodel"); 215 218 } 219 220 216 221 /*i/o: */ 217 222 IoModelFetchData(&this->io_gather,iomodel_handle,"io_gather"); … … 242 247 this->z=NULL; 243 248 this->elements=NULL; 249 this->elementconnectivity=NULL; 244 250 this->elements_type=NULL; 245 251 this->vertices_type=NULL; … … 275 281 this->surface=NULL; 276 282 this->bed=NULL; 283 this->bathymetry=NULL; 277 284 this->elementoniceshelf=NULL; 278 285 this->elementonwater=NULL; -
issm/trunk/src/c/objects/IoModel.h
r6946 r7089 32 32 double* elements_type; 33 33 double* vertices_type; 34 double* elementconnectivity; 34 35 35 36 /*3d mesh: */ … … 83 84 double* surface; 84 85 double* bed; 86 double* bathymetry; 85 87 double* elementoniceshelf; 86 88 double* elementonwater; … … 189 191 int qmu_save_femmodel; 190 192 193 /*grounding line migration: */ 194 int grounding_line_migration; 195 191 196 /*exterior partitioning data, to be carried around: */ 192 197 bool* my_elements; -
issm/trunk/src/c/objects/Loads/Riftfront.cpp
r6748 r7089 87 87 this->length=*(iomodel->riftinfo+RIFTINFOSIZE*i+6); 88 88 this->fraction=*(iomodel->riftinfo+RIFTINFOSIZE*i+9); 89 this->state= *(iomodel->riftinfo+RIFTINFOSIZE*i+11);89 this->state=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+11); 90 90 91 91 //intialize inputs, and add as many inputs per element as requested: -
issm/trunk/src/c/objects/Node.h
r6412 r7089 21 21 class Node: public Object ,public Update{ 22 22 23 p rivate:23 public: 24 24 25 25 int id; //unique arbitrary id. … … 31 31 int analysis_type; 32 32 33 public:34 33 35 34 /*Node constructors, destructors {{{1*/ -
issm/trunk/src/c/objects/objects.h
r6200 r7089 57 57 #include "./ElementResults/DoubleElementResult.h" 58 58 #include "./ElementResults/TriaVertexElementResult.h" 59 #include "./ElementResults/PentaVertexElementResult.h" 59 #include "./ElementResults/PentaVertexElementResult.h" 60 #include "./ElementResults/BoolElementResult.h" 60 61 61 62 /*ExternalResults: */ -
issm/trunk/src/c/solutions/transient2d_core.cpp
r6962 r7089 22 22 bool time_adapt=false; 23 23 int output_frequency; 24 bool grounding_line_migration; 24 25 25 26 /*intermediary: */ … … 35 36 femmodel->parameters->FindParam(&output_frequency,OutputFrequencyEnum); 36 37 femmodel->parameters->FindParam(&time_adapt,TimeAdaptEnum); 38 femmodel->parameters->FindParam(&grounding_line_migration,GroundingLineMigrationEnum); 37 39 38 40 /*initialize: */ … … 58 60 _printf_(VerboseSolution(),"%s\n"," computing new thickness"); 59 61 prognostic_core(femmodel); 60 62 63 if (grounding_line_migration){ 64 _printf_(VerboseSolution(),"%s\n"," computing new grounding line position"); 65 GroundingLineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); 66 } 67 61 68 if(solution_type==Transient2DSolutionEnum && !control_analysis && (step%output_frequency==0 || time==finaltime)){ 62 69 _printf_(VerboseSolution(),"%s\n"," saving results\n"); … … 68 75 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); 69 76 InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); 77 if(grounding_line_migration)InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ElementOnIceShelfEnum,step,time); 70 78 71 79 /*unload results*/
Note:
See TracChangeset
for help on using the changeset viewer.