![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 3422
- Timestamp:
- 04/07/10 15:30:57 (15 years ago)
- Location:
- issm/trunk/src/c
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp
r3417 r3422 15 15 void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ 16 16 17 18 /*output: int* epart, int* my_grids, double* my_bordergrids*/ 19 20 17 /*Intermediary*/ 21 18 int i,j,k,n; 22 extern int my_rank;23 extern int num_procs;24 19 25 20 /*DataSets: */ … … 35 30 Matice* matice = NULL; 36 31 Matpar* matpar = NULL; 37 ElementProperties* penta_properties=NULL;38 32 39 /*output: */ 40 int* epart=NULL; //element partitioning. 41 int* npart=NULL; //node partitioning. 42 int* my_grids=NULL; 43 double* my_bordergrids=NULL; 44 45 46 /*intermediary: */ 47 int elements_width; //size of elements 48 double B_avg; 49 50 /*matice constructor input: */ 51 int matice_mid; 52 double matice_B; 53 double matice_n; 54 55 /*penta constructor input: */ 56 57 int penta_id; 58 int penta_matice_id; 59 int penta_matpar_id; 60 int penta_numpar_id; 61 int penta_node_ids[6]; 62 double penta_h[6]; 63 double penta_s[6]; 64 double penta_b[6]; 65 double penta_k[6]; 66 int penta_friction_type; 67 double penta_p; 68 double penta_q; 69 int penta_shelf; 70 int penta_onbed; 71 int penta_onsurface; 72 int penta_collapse=0; 73 double penta_melting[6]; 74 double penta_accumulation[6]; 75 bool penta_onwater; 76 77 /*matpar constructor input: */ 78 int matpar_mid; 79 double matpar_rho_ice; 80 double matpar_rho_water; 81 double matpar_heatcapacity; 82 double matpar_thermalconductivity; 83 double matpar_latentheat; 84 double matpar_beta; 85 double matpar_meltingpoint; 86 double matpar_mixed_layer_capacity; 87 double matpar_thermal_exchange_velocity; 88 double matpar_g; 89 90 /* node constructor input: */ 91 int node_id; 92 int vertex_id; 93 int node_partitionborder=0; 94 int node_onbed; 95 int node_onsurface; 96 int node_onshelf; 97 int node_onsheet; 98 int node_upper_node_id; 99 int node_numdofs; 100 101 102 #ifdef _PARALLEL_ 103 /*Metis partitioning: */ 104 int range; 105 Vec gridborder=NULL; 106 int my_numgrids; 107 int* all_numgrids=NULL; 108 int gridcount; 109 int count; 110 #endif 111 int first_grid_index; 33 /*Now, do we have Stokes elements?*/ 34 if (strcmp(iomodel->meshtype,"2d")==0) ISSMERROR("Stokes elements only supported in 3d!"); 35 if (!iomodel->isstokes)goto cleanup_and_return; 112 36 113 37 /*First create the elements, nodes and material properties: */ … … 117 41 materials = new DataSet(MaterialsEnum()); 118 42 119 /* Now, is the flag ishutter on? otherwise, do nothing: */120 if (!iomodel->isstokes)goto cleanup_and_return;43 /*Partition elements and vertices and nodes: */ 44 Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); 121 45 122 /*Width of elements: */ 123 if(strcmp(iomodel->meshtype,"2d")==0){ 124 elements_width=3; //tria elements 125 } 126 else{ 127 elements_width=6; //penta elements 128 } 46 /*Fetch data needed: */ 47 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 48 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); 49 IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); 50 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); 51 IoModelFetchData(&iomodel->drag,NULL,NULL,iomodel_handle,"drag"); 52 IoModelFetchData(&iomodel->p,NULL,NULL,iomodel_handle,"p"); 53 IoModelFetchData(&iomodel->q,NULL,NULL,iomodel_handle,"q"); 54 IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); 55 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); 56 IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); 57 IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); 58 IoModelFetchData(&iomodel->B,NULL,NULL,iomodel_handle,"B"); 59 IoModelFetchData(&iomodel->n,NULL,NULL,iomodel_handle,"n"); 60 IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation"); 61 IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); 62 IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); 129 63 64 for (i=0;i<iomodel->numberofelements;i++){ 130 65 131 #ifdef _PARALLEL_ 132 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ 133 if(strcmp(iomodel->meshtype,"2d")==0){ 134 /*load elements: */ 135 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); 136 } 137 else{ 138 /*load elements2d: */ 139 IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); 140 } 66 if(iomodel->my_elements[i]){ 141 67 68 /*Create and add penta element to elements dataset: */ 69 elements->AddObject(new Penta(i,iomodel)); 142 70 143 MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); 71 /*Create and add material property to materials dataset: */ 72 materials->AddObject(new Matice(i,iomodel,6)); 144 73 145 /*Free elements and elements2d: */146 xfree((void**)&iomodel->elements);147 xfree((void**)&iomodel->elements2d);148 149 /*Used later on: */150 my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int));151 #endif152 153 154 /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */155 156 /*2d mesh: */157 if (strcmp(iomodel->meshtype,"2d")==0){158 ISSMERROR(" stokes elements only supported in 3d!");159 }160 else{ // if (strcmp(meshtype,"2d")==0)161 162 /*Fetch data needed: */163 IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements");164 IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness");165 IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface");166 IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed");167 IoModelFetchData(&iomodel->drag,NULL,NULL,iomodel_handle,"drag");168 IoModelFetchData(&iomodel->p,NULL,NULL,iomodel_handle,"p");169 IoModelFetchData(&iomodel->q,NULL,NULL,iomodel_handle,"q");170 IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf");171 IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed");172 IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface");173 IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type");174 IoModelFetchData(&iomodel->B,NULL,NULL,iomodel_handle,"B");175 IoModelFetchData(&iomodel->n,NULL,NULL,iomodel_handle,"n");176 IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation");177 IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting");178 IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater");179 180 for (i=0;i<iomodel->numberofelements;i++){181 #ifdef _PARALLEL_182 /*We are using our element partition to decide which elements will be created on this node: */183 if(my_rank==epart[i]){184 #endif185 186 187 /*name and id: */188 penta_id=i+1; //matlab indexing.189 penta_matice_id=i+1; //refers to the corresponding material property card190 penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card191 penta_numpar_id=1;192 193 /*vertices,thickness,surface,bed and drag: */194 for(j=0;j<6;j++){195 penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j);196 penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1));197 penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1));198 penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1));199 penta_k[j]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+j)-1));200 penta_melting[j]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+j)-1));201 penta_accumulation[j]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+j)-1));202 }203 204 /*basal drag:*/205 penta_friction_type=(int)iomodel->drag_type;206 207 penta_p=iomodel->p[i];208 penta_q=iomodel->q[i];209 210 /*diverse: */211 penta_shelf=(int)*(iomodel->elementoniceshelf+i);212 penta_onbed=(int)*(iomodel->elementonbed+i);213 penta_onsurface=(int)*(iomodel->elementonsurface+i);214 penta_onwater=(bool)*(iomodel->elementonwater+i);215 penta_collapse=0;216 217 if (*(iomodel->elements_type+2*i+1)==StokesFormulationEnum()){218 219 /*Create properties: */220 penta_properties=new ElementProperties(6,penta_h, penta_s, penta_b, penta_k, penta_melting, penta_accumulation, NULL, penta_friction_type, penta_p, penta_q, penta_shelf, penta_onbed, penta_onwater, penta_onsurface, penta_collapse, UNDEF);221 222 /*Create Penta using its constructor:*/223 penta=new Penta(penta_id,penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties);224 225 /*Delete properties: */226 delete penta_properties;227 228 /*Add penta element to elements dataset: */229 elements->AddObject(penta);230 }231 232 233 /*Deal with material:*/234 matice_mid=i+1; //same as the material id from the geom2 elements.235 /*Average B over 6 element grids: */236 B_avg=0;237 for(j=0;j<6;j++){238 B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1));239 }240 B_avg=B_avg/6;241 matice_B= B_avg;242 matice_n=(double)*(iomodel->n+i);243 244 if (*(iomodel->elements_type+2*i+1)==StokesFormulationEnum()){245 246 /*Create matice using its constructor:*/247 matice= new Matice(matice_mid,matice_B,matice_n);248 249 /*Add matice element to materials dataset: */250 materials->AddObject(matice);251 }252 253 #ifdef _PARALLEL_254 /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use255 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)256 into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids257 will hold which grids belong to this partition*/258 my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1;259 my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1;260 my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1;261 my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1;262 my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1;263 my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1;264 #endif265 266 #ifdef _PARALLEL_267 }//if(my_rank==epart[i])268 #endif269 270 }//for (i=0;i<numberofelements;i++)271 272 /*Free data: */273 xfree((void**)&iomodel->elements);274 xfree((void**)&iomodel->thickness);275 xfree((void**)&iomodel->surface);276 xfree((void**)&iomodel->bed);277 xfree((void**)&iomodel->drag);278 xfree((void**)&iomodel->p);279 xfree((void**)&iomodel->q);280 xfree((void**)&iomodel->elementoniceshelf);281 xfree((void**)&iomodel->elementonbed);282 xfree((void**)&iomodel->elementonsurface);283 xfree((void**)&iomodel->elements_type);284 xfree((void**)&iomodel->n);285 xfree((void**)&iomodel->B);286 xfree((void**)&iomodel->accumulation);287 xfree((void**)&iomodel->melting);288 xfree((void**)&iomodel->elementonwater);289 290 } //if (strcmp(meshtype,"2d")==0)291 292 /*Add one constant material property to materials: */293 matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials294 matpar_g=iomodel->g;295 matpar_rho_ice=iomodel->rho_ice;296 matpar_rho_water=iomodel->rho_water;297 matpar_thermalconductivity=iomodel->thermalconductivity;298 matpar_heatcapacity=iomodel->heatcapacity;299 matpar_latentheat=iomodel->latentheat;300 matpar_beta=iomodel->beta;301 matpar_meltingpoint=iomodel->meltingpoint;302 matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;303 matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;304 305 /*Create matpar object using its constructor: */306 matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity,307 matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity,308 matpar_thermal_exchange_velocity,matpar_g);309 310 /*Add to materials datset: */311 materials->AddObject(matpar);312 313 #ifdef _PARALLEL_314 /*From the element partitioning, we can determine which grids are on the inside of this cpu's315 *element partition, and which are on its border with other nodes:*/316 gridborder=NewVec(iomodel->numberofnodes);317 318 for (i=0;i<iomodel->numberofnodes;i++){319 if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES);320 }321 VecAssemblyBegin(gridborder);322 VecAssemblyEnd(gridborder);323 324 #ifdef _ISSM_DEBUG_325 VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD);326 #endif327 328 VecToMPISerial(&my_bordergrids,gridborder);329 330 #ifdef _ISSM_DEBUG_331 if(my_rank==0){332 for (i=0;i<iomodel->numberofnodes;i++){333 printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]);334 }335 }336 #endif337 #endif338 339 /*Partition penalties in 3d: */340 if(strcmp(iomodel->meshtype,"3d")==0){341 342 /*Get penalties: */343 IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties");344 345 if(iomodel->numpenalties){346 347 iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int));348 #ifdef _SERIAL_349 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=1;350 #else351 for(i=0;i<iomodel->numpenalties;i++)iomodel->penaltypartitioning[i]=-1;352 353 for(i=0;i<iomodel->numpenalties;i++){354 first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1);355 if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids356 /*All grids that are being penalised belong to this node's internal grid partition.:*/357 iomodel->penaltypartitioning[i]=1;358 }359 if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border360 iomodel->penaltypartitioning[i]=0;361 }362 }363 #endif364 74 } 365 75 366 /*Free penalties: */ 367 xfree((void**)&iomodel->penalties); 368 } 76 }//for (i=0;i<numberofelements;i++) 369 77 370 /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. 371 We can therefore determine which grids are internal to this node's partition 372 and which ones are shared with other nodes because they are on the border of this node's partition. Knowing 373 that, go and create the grids*/ 78 /*Free data: */ 79 xfree((void**)&iomodel->elements); 80 xfree((void**)&iomodel->thickness); 81 xfree((void**)&iomodel->surface); 82 xfree((void**)&iomodel->bed); 83 xfree((void**)&iomodel->drag); 84 xfree((void**)&iomodel->p); 85 xfree((void**)&iomodel->q); 86 xfree((void**)&iomodel->elementoniceshelf); 87 xfree((void**)&iomodel->elementonbed); 88 xfree((void**)&iomodel->elementonsurface); 89 xfree((void**)&iomodel->elements_type); 90 xfree((void**)&iomodel->n); 91 xfree((void**)&iomodel->B); 92 xfree((void**)&iomodel->accumulation); 93 xfree((void**)&iomodel->melting); 94 xfree((void**)&iomodel->elementonwater); 374 95 375 /*Create nodes from x,y,z, as well as the spc values on those grids: */ 376 96 /*Add new constrant material property to materials, at the end: */ 97 materials->AddObject(new Matpar(iomodel)); 98 377 99 /*First fetch data: */ 378 if (strcmp(iomodel->meshtype,"3d")==0){ 379 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 380 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 381 } 100 IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); 101 IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); 382 102 IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); 383 103 IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); … … 392 112 IoModelFetchData(&iomodel->borderstokes,NULL,NULL,iomodel_handle,"borderstokes"); 393 113 114 for (i=0;i<iomodel->numberofvertices;i++){ 394 115 395 /*Get number of dofs per node: */396 DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type);116 /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ 117 if(iomodel->my_vertices[i]){ 397 118 398 for (i=0;i<iomodel->numberofnodes;i++){ 399 #ifdef _PARALLEL_ 400 /*keep only this partition's nodes:*/ 401 if((my_grids[i]==1)){ 402 #endif 119 /*Add vertex to vertices dataset: */ 120 vertices->AddObject(new Vertex(i,iomodel)); 403 121 404 /*create vertex: */ 405 vertex_id=i+1; 406 vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); 122 /*Add node to nodes dataset: */ 123 nodes->AddObject(new Node(i,iomodel)); 407 124 408 vertices->AddObject(vertex);409 410 node_id=i+1; //matlab indexing411 412 #ifdef _PARALLEL_413 if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border414 node_partitionborder=1;415 125 } 416 else{417 node_partitionborder=0;418 }419 #else420 node_partitionborder=0;421 #endif422 423 node_properties.SetProperties(424 (int)iomodel->gridonbed[i],425 (int)iomodel->gridonsurface[i],426 (int)iomodel->gridoniceshelf[i],427 (int)iomodel->gridonicesheet[i]);428 429 430 if (strcmp(iomodel->meshtype,"3d")==0){431 if (isnan(iomodel->uppernodes[i])){432 node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves.433 }434 else{435 node_upper_node_id=(int)iomodel->uppernodes[i];436 }437 }438 else{439 /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/440 node_upper_node_id=node_id;441 }442 443 /*Create node using its constructor: */444 node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties);445 446 /*set single point constraints.: */447 /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */448 if (iomodel->borderstokes[i]){449 //freeze everything except pressure450 node->FreezeDof(1);451 node->FreezeDof(2);452 node->FreezeDof(3);453 }454 else if (iomodel->gridonstokes[i]==0){455 for(k=1;k<=node_numdofs;k++){456 node->FreezeDof(k);457 }458 }459 460 461 /*Add node to nodes dataset: */462 nodes->AddObject(node);463 464 #ifdef _PARALLEL_465 } //if((my_grids[i]==1))466 #endif467 126 } 468 469 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these470 * datasets, it will not be redone: */471 elements->Presort();472 nodes->Presort();473 vertices->Presort();474 materials->Presort();475 127 476 128 /*Clean fetched data: */ … … 489 141 xfree((void**)&iomodel->gridoniceshelf); 490 142 491 492 /*Keep partitioning information into iomodel*/ 493 iomodel->epart=epart; 494 iomodel->my_grids=my_grids; 495 iomodel->my_bordergrids=my_bordergrids; 496 497 /*Free ressources:*/ 498 #ifdef _PARALLEL_ 499 xfree((void**)&all_numgrids); 500 xfree((void**)&npart); 501 VecFree(&gridborder); 502 #endif 143 /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these 144 * datasets, it will not be redone: */ 145 elements->Presort(); 146 nodes->Presort(); 147 vertices->Presort(); 148 materials->Presort(); 503 149 504 150 cleanup_and_return: -
issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp
r3421 r3422 83 83 xfree((void**)&iomodel->elementonwater); 84 84 85 /*Add new constrant material property t go materials, at the end: */85 /*Add new constrant material property to materials, at the end: */ 86 86 materials->AddObject(new Matpar(iomodel)); 87 87 -
issm/trunk/src/c/objects/Node.cpp
r3420 r3422 62 62 int upper_node_id; 63 63 64 65 64 /*id: */ 66 65 this->id=i+1; //matlab indexing … … 97 96 this->hvertex.Init(&vertex_id,1); //node id is the same as the vertex id, continuous galerkin! 98 97 this->hupper_node.Init(&upper_node_id,1); 99 100 98 101 99 /*set single point constraints.: */ … … 109 107 } 110 108 if (iomodel->gridonhutter[i]){ 109 for(k=1;k<=numdofs;k++){ 110 this->FreezeDof(k); 111 } 112 } 113 /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ 114 if (iomodel->borderstokes[i]){ 115 //freeze everything except pressure 116 node->FreezeDof(1); 117 node->FreezeDof(2); 118 node->FreezeDof(3); 119 } 120 else if (iomodel->gridonstokes[i]==0){ 111 121 for(k=1;k<=numdofs;k++){ 112 122 this->FreezeDof(k);
Note:
See TracChangeset
for help on using the changeset viewer.
