Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 3421) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 3422) @@ -15,11 +15,6 @@ void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - - /*output: int* epart, int* my_grids, double* my_bordergrids*/ - - + /*Intermediary*/ int i,j,k,n; - extern int my_rank; - extern int num_procs; /*DataSets: */ @@ -35,79 +30,8 @@ Matice* matice = NULL; Matpar* matpar = NULL; - ElementProperties* penta_properties=NULL; - /*output: */ - int* epart=NULL; //element partitioning. - int* npart=NULL; //node partitioning. - int* my_grids=NULL; - double* my_bordergrids=NULL; - - - /*intermediary: */ - int elements_width; //size of elements - double B_avg; - - /*matice constructor input: */ - int matice_mid; - double matice_B; - double matice_n; - - /*penta constructor input: */ - - int penta_id; - int penta_matice_id; - int penta_matpar_id; - int penta_numpar_id; - int penta_node_ids[6]; - double penta_h[6]; - double penta_s[6]; - double penta_b[6]; - double penta_k[6]; - int penta_friction_type; - double penta_p; - double penta_q; - int penta_shelf; - int penta_onbed; - int penta_onsurface; - int penta_collapse=0; - double penta_melting[6]; - double penta_accumulation[6]; - bool penta_onwater; - - /*matpar constructor input: */ - int matpar_mid; - double matpar_rho_ice; - double matpar_rho_water; - double matpar_heatcapacity; - double matpar_thermalconductivity; - double matpar_latentheat; - double matpar_beta; - double matpar_meltingpoint; - double matpar_mixed_layer_capacity; - double matpar_thermal_exchange_velocity; - double matpar_g; - - /* node constructor input: */ - int node_id; - int vertex_id; - int node_partitionborder=0; - int node_onbed; - int node_onsurface; - int node_onshelf; - int node_onsheet; - int node_upper_node_id; - int node_numdofs; - - - #ifdef _PARALLEL_ - /*Metis partitioning: */ - int range; - Vec gridborder=NULL; - int my_numgrids; - int* all_numgrids=NULL; - int gridcount; - int count; - #endif - int first_grid_index; + /*Now, do we have Stokes elements?*/ + if (strcmp(iomodel->meshtype,"2d")==0) ISSMERROR("Stokes elements only supported in 3d!"); + if (!iomodel->isstokes)goto cleanup_and_return; /*First create the elements, nodes and material properties: */ @@ -117,267 +41,63 @@ materials = new DataSet(MaterialsEnum()); - /*Now, is the flag ishutter on? otherwise, do nothing: */ - if (!iomodel->isstokes)goto cleanup_and_return; + /*Partition elements and vertices and nodes: */ + Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); - /*Width of elements: */ - if(strcmp(iomodel->meshtype,"2d")==0){ - elements_width=3; //tria elements - } - else{ - elements_width=6; //penta elements - } + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->drag,NULL,NULL,iomodel_handle,"drag"); + IoModelFetchData(&iomodel->p,NULL,NULL,iomodel_handle,"p"); + IoModelFetchData(&iomodel->q,NULL,NULL,iomodel_handle,"q"); + IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); + IoModelFetchData(&iomodel->B,NULL,NULL,iomodel_handle,"B"); + IoModelFetchData(&iomodel->n,NULL,NULL,iomodel_handle,"n"); + IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation"); + IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + for (i=0;inumberofelements;i++){ - #ifdef _PARALLEL_ - /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ - if(strcmp(iomodel->meshtype,"2d")==0){ - /*load elements: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - } - else{ - /*load elements2d: */ - IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); - } + if(iomodel->my_elements[i]){ + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Penta(i,iomodel)); - MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,6)); - /*Free elements and elements2d: */ - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->elements2d); - - /*Used later on: */ - my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); - #endif - - - /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ - - /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ - ISSMERROR(" stokes elements only supported in 3d!"); - } - else{ // if (strcmp(meshtype,"2d")==0) - - /*Fetch data needed: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); - IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); - IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); - IoModelFetchData(&iomodel->drag,NULL,NULL,iomodel_handle,"drag"); - IoModelFetchData(&iomodel->p,NULL,NULL,iomodel_handle,"p"); - IoModelFetchData(&iomodel->q,NULL,NULL,iomodel_handle,"q"); - IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); - IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); - IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); - IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); - IoModelFetchData(&iomodel->B,NULL,NULL,iomodel_handle,"B"); - IoModelFetchData(&iomodel->n,NULL,NULL,iomodel_handle,"n"); - IoModelFetchData(&iomodel->accumulation,NULL,NULL,iomodel_handle,"accumulation"); - IoModelFetchData(&iomodel->melting,NULL,NULL,iomodel_handle,"melting"); - IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); - - for (i=0;inumberofelements;i++){ - #ifdef _PARALLEL_ - /*We are using our element partition to decide which elements will be created on this node: */ - if(my_rank==epart[i]){ - #endif - - - /*name and id: */ - penta_id=i+1; //matlab indexing. - penta_matice_id=i+1; //refers to the corresponding material property card - penta_matpar_id=iomodel->numberofelements+1;//refers to the corresponding parmat property card - penta_numpar_id=1; - - /*vertices,thickness,surface,bed and drag: */ - for(j=0;j<6;j++){ - penta_node_ids[j]=(int)*(iomodel->elements+elements_width*i+j); - penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_k[j]=*(iomodel->drag+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_melting[j]=*(iomodel->melting+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - penta_accumulation[j]=*(iomodel->accumulation+ ((int)*(iomodel->elements+elements_width*i+j)-1)); - } - - /*basal drag:*/ - penta_friction_type=(int)iomodel->drag_type; - - penta_p=iomodel->p[i]; - penta_q=iomodel->q[i]; - - /*diverse: */ - penta_shelf=(int)*(iomodel->elementoniceshelf+i); - penta_onbed=(int)*(iomodel->elementonbed+i); - penta_onsurface=(int)*(iomodel->elementonsurface+i); - penta_onwater=(bool)*(iomodel->elementonwater+i); - penta_collapse=0; - - if (*(iomodel->elements_type+2*i+1)==StokesFormulationEnum()){ - - /*Create properties: */ - 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); - - /*Create Penta using its constructor:*/ - penta=new Penta(penta_id,penta_node_ids, penta_matice_id, penta_matpar_id, penta_numpar_id, penta_properties); - - /*Delete properties: */ - delete penta_properties; - - /*Add penta element to elements dataset: */ - elements->AddObject(penta); - } - - - /*Deal with material:*/ - matice_mid=i+1; //same as the material id from the geom2 elements. - /*Average B over 6 element grids: */ - B_avg=0; - for(j=0;j<6;j++){ - B_avg+=*(iomodel->B+((int)*(iomodel->elements+elements_width*i+j)-1)); - } - B_avg=B_avg/6; - matice_B= B_avg; - matice_n=(double)*(iomodel->n+i); - - if (*(iomodel->elements_type+2*i+1)==StokesFormulationEnum()){ - - /*Create matice using its constructor:*/ - matice= new Matice(matice_mid,matice_B,matice_n); - - /*Add matice element to materials dataset: */ - materials->AddObject(matice); - } - - #ifdef _PARALLEL_ - /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use - *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) - into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids - will hold which grids belong to this partition*/ - my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1; - my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1; - #endif - - #ifdef _PARALLEL_ - }//if(my_rank==epart[i]) - #endif - - }//for (i=0;ielements); - xfree((void**)&iomodel->thickness); - xfree((void**)&iomodel->surface); - xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->drag); - xfree((void**)&iomodel->p); - xfree((void**)&iomodel->q); - xfree((void**)&iomodel->elementoniceshelf); - xfree((void**)&iomodel->elementonbed); - xfree((void**)&iomodel->elementonsurface); - xfree((void**)&iomodel->elements_type); - xfree((void**)&iomodel->n); - xfree((void**)&iomodel->B); - xfree((void**)&iomodel->accumulation); - xfree((void**)&iomodel->melting); - xfree((void**)&iomodel->elementonwater); - - } //if (strcmp(meshtype,"2d")==0) - - /*Add one constant material property to materials: */ - matpar_mid=iomodel->numberofelements+1; //put it at the end of the materials - matpar_g=iomodel->g; - matpar_rho_ice=iomodel->rho_ice; - matpar_rho_water=iomodel->rho_water; - matpar_thermalconductivity=iomodel->thermalconductivity; - matpar_heatcapacity=iomodel->heatcapacity; - matpar_latentheat=iomodel->latentheat; - matpar_beta=iomodel->beta; - matpar_meltingpoint=iomodel->meltingpoint; - matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; - matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; - - /*Create matpar object using its constructor: */ - matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, - matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, - matpar_thermal_exchange_velocity,matpar_g); - - /*Add to materials datset: */ - materials->AddObject(matpar); - - #ifdef _PARALLEL_ - /*From the element partitioning, we can determine which grids are on the inside of this cpu's - *element partition, and which are on its border with other nodes:*/ - gridborder=NewVec(iomodel->numberofnodes); - - for (i=0;inumberofnodes;i++){ - if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); - } - VecAssemblyBegin(gridborder); - VecAssemblyEnd(gridborder); - - #ifdef _ISSM_DEBUG_ - VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); - #endif - - VecToMPISerial(&my_bordergrids,gridborder); - - #ifdef _ISSM_DEBUG_ - if(my_rank==0){ - for (i=0;inumberofnodes;i++){ - printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); - } - } - #endif - #endif - - /*Partition penalties in 3d: */ - if(strcmp(iomodel->meshtype,"3d")==0){ - - /*Get penalties: */ - IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); - - if(iomodel->numpenalties){ - - iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int)); - #ifdef _SERIAL_ - for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=1; - #else - for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=-1; - - for(i=0;inumpenalties;i++){ - first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1); - if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids - /*All grids that are being penalised belong to this node's internal grid partition.:*/ - iomodel->penaltypartitioning[i]=1; - } - if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border - iomodel->penaltypartitioning[i]=0; - } - } - #endif } - /*Free penalties: */ - xfree((void**)&iomodel->penalties); - } + }//for (i=0;ielements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->drag); + xfree((void**)&iomodel->p); + xfree((void**)&iomodel->q); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elements_type); + xfree((void**)&iomodel->n); + xfree((void**)&iomodel->B); + xfree((void**)&iomodel->accumulation); + xfree((void**)&iomodel->melting); + xfree((void**)&iomodel->elementonwater); - /*Create nodes from x,y,z, as well as the spc values on those grids: */ - + /*Add new constrant material property to materials, at the end: */ + materials->AddObject(new Matpar(iomodel)); + /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ - IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); - IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); - } + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); @@ -392,85 +112,17 @@ IoModelFetchData(&iomodel->borderstokes,NULL,NULL,iomodel_handle,"borderstokes"); + for (i=0;inumberofvertices;i++){ - /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); + /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ + if(iomodel->my_vertices[i]){ - for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if((my_grids[i]==1)){ - #endif + /*Add vertex to vertices dataset: */ + vertices->AddObject(new Vertex(i,iomodel)); - /*create vertex: */ - vertex_id=i+1; - vertex=new Vertex(vertex_id, iomodel->x[i],iomodel->y[i],iomodel->z[i],(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i])); + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(i,iomodel)); - vertices->AddObject(vertex); - - node_id=i+1; //matlab indexing - - #ifdef _PARALLEL_ - if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border - node_partitionborder=1; } - else{ - node_partitionborder=0; - } - #else - node_partitionborder=0; - #endif - - node_properties.SetProperties( - (int)iomodel->gridonbed[i], - (int)iomodel->gridonsurface[i], - (int)iomodel->gridoniceshelf[i], - (int)iomodel->gridonicesheet[i]); - - - if (strcmp(iomodel->meshtype,"3d")==0){ - if (isnan(iomodel->uppernodes[i])){ - node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. - } - else{ - node_upper_node_id=(int)iomodel->uppernodes[i]; - } - } - else{ - /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ - node_upper_node_id=node_id; - } - - /*Create node using its constructor: */ - node=new Node(node_id,vertex_id, node_upper_node_id, node_partitionborder,node_numdofs, &node_properties); - - /*set single point constraints.: */ - /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ - if (iomodel->borderstokes[i]){ - //freeze everything except pressure - node->FreezeDof(1); - node->FreezeDof(2); - node->FreezeDof(3); - } - else if (iomodel->gridonstokes[i]==0){ - for(k=1;k<=node_numdofs;k++){ - node->FreezeDof(k); - } - } - - - /*Add node to nodes dataset: */ - nodes->AddObject(node); - - #ifdef _PARALLEL_ - } //if((my_grids[i]==1)) - #endif } - - /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these - * datasets, it will not be redone: */ - elements->Presort(); - nodes->Presort(); - vertices->Presort(); - materials->Presort(); /*Clean fetched data: */ @@ -489,16 +141,10 @@ xfree((void**)&iomodel->gridoniceshelf); - - /*Keep partitioning information into iomodel*/ - iomodel->epart=epart; - iomodel->my_grids=my_grids; - iomodel->my_bordergrids=my_bordergrids; - - /*Free ressources:*/ - #ifdef _PARALLEL_ - xfree((void**)&all_numgrids); - xfree((void**)&npart); - VecFree(&gridborder); - #endif + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + elements->Presort(); + nodes->Presort(); + vertices->Presort(); + materials->Presort(); cleanup_and_return: Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3421) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3422) @@ -83,5 +83,5 @@ xfree((void**)&iomodel->elementonwater); - /*Add new constrant material property tgo materials, at the end: */ + /*Add new constrant material property to materials, at the end: */ materials->AddObject(new Matpar(iomodel)); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3421) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3422) @@ -62,5 +62,4 @@ int upper_node_id; - /*id: */ this->id=i+1; //matlab indexing @@ -97,5 +96,4 @@ this->hvertex.Init(&vertex_id,1); //node id is the same as the vertex id, continuous galerkin! this->hupper_node.Init(&upper_node_id,1); - /*set single point constraints.: */ @@ -109,4 +107,16 @@ } if (iomodel->gridonhutter[i]){ + for(k=1;k<=numdofs;k++){ + this->FreezeDof(k); + } + } + /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ + if (iomodel->borderstokes[i]){ + //freeze everything except pressure + node->FreezeDof(1); + node->FreezeDof(2); + node->FreezeDof(3); + } + else if (iomodel->gridonstokes[i]==0){ for(k=1;k<=numdofs;k++){ this->FreezeDof(k);