Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3422) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3423) @@ -23,12 +23,4 @@ DataSet* materials = NULL; - /*Objects: */ - Node* node = NULL; - Vertex* vertex = NULL; - Tria* tria = NULL; - Penta* penta = NULL; - Matice* matice = NULL; - Matpar* matpar = NULL; - /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ if (!iomodel->ismacayealpattyn)goto cleanup_and_return; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 3422) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 3423) @@ -2,5 +2,4 @@ * CreateElementsNodesAndMaterialsDiagnosticHutter.c: */ - #include "../../DataSet/DataSet.h" @@ -13,187 +12,17 @@ #include "../IoModel.h" +void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ -void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - - - int i,k; - extern int my_rank; - extern int num_procs; + /*Intermediary*/ + int i,j,k,n; /*DataSets: */ DataSet* elements = NULL; DataSet* nodes = NULL; + DataSet* vertices = NULL; DataSet* materials = NULL; - - /*Objects: */ - Node* node = NULL; - Matice* matice = NULL; - Matpar* matpar = NULL; - Beam* beam = NULL; - Sing* sing = NULL; - ElementProperties* sing_properties=NULL; - ElementProperties* beam_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; - - /*sing constructor input: */ - int sing_id; - int sing_matice_id; - int sing_matpar_id; - int sing_numpar_id; - int sing_node_ids[1]; - double sing_h[1],sing_k[1]; - - /*beam constructor input: */ - int beam_id; - int beam_matice_id; - int beam_matpar_id; - int beam_numpar_id; - int beam_node_ids[2]; - double beam_h[2]; - double beam_s[2]; - double beam_b[2]; - double beam_k[2]; - bool beam_onbed; - - /*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 node_partitionborder=0; - double node_x[3]; - double node_sigma; - 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 - - /*First create the elements, nodes and material properties: */ - elements = new DataSet(ElementsEnum()); - nodes = new DataSet(NodesEnum()); - materials = new DataSet(MaterialsEnum()); /*Now, is the flag ishutter on? otherwise, do nothing: */ if (!iomodel->ishutter)goto cleanup_and_return; - - /*Width of elements: */ - if(strcmp(iomodel->meshtype,"2d")==0){ - elements_width=3; //sing elements - } - else{ - elements_width=6; //beam elements - } - - #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"); - } - - MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); - - /*Free elements and elements2d: */ - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->elements2d); - - /*Used later on: */ - my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); - #endif - - #ifdef _PARALLEL_ - if(strcmp(iomodel->meshtype,"2d")==0){ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - for (i=0;inumberofelements;i++){ - if(my_rank==epart[i]){ - /*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; - } - } - } - else{ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - for (i=0;inumberofelements;i++){ - if(my_rank==epart[i]){ - 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; - } - } - } - - /*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 /*Hutter elements can be partitioned using epart, even if @@ -216,47 +45,14 @@ for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if(my_grids[i]){ - #endif - if(iomodel->gridonhutter[i]){ - - /*Deal with sing element: */ - sing_id=i+1; - sing_matice_id=i+1; //refers to the corresponding material property card - sing_matpar_id=iomodel->numberofnodes+1;//refers to the corresponding matpar property card - sing_numpar_id=1; - sing_node_ids[0]=i+1; - sing_h[0]=iomodel->thickness[i]; - sing_k[0]=iomodel->drag[i]; + if(iomodel->my_nodes[i]){ - /*Create properties: */ - sing_properties=new ElementProperties(1,sing_h,NULL,NULL, sing_k,NULL,NULL, NULL, UNDEF, (double)UNDEF, (double)UNDEF, UNDEF, UNDEF,true, UNDEF,UNDEF,UNDEF); + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Sing(i,iomodel)); - /*Create sing element using its constructor:*/ - sing=new Sing(sing_id, sing_node_ids, sing_matice_id, sing_matpar_id, sing_numpar_id, sing_properties); + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,1)); - /*delete properties: */ - delete sing_properties; - - /*Add sing element to elements dataset: */ - elements->AddObject(sing); - - /*Deal with material property card: */ - matice_mid=i+1; //same as the material id from the geom2 elements. - matice_B=iomodel->B[i]; - matice_n=(double)iomodel->n[1]; //n defined on elements not grids, so take the first value everywhere - - /*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_ - } //if(my_rank==npart[i]) - #endif } //for (i=0;inumberofnodes;i++) @@ -265,59 +61,18 @@ for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if(my_grids[i]){ - #endif - if(iomodel->gridonhutter[i]){ + if(iomodel->my_nodes[i]){ + if(iomodel->gridonhutter[i]){ + if(!iomodel->gridonsurface[i]){ - if(!iomodel->gridonsurface[i]){ - - /*Deal with sing element: */ - beam_id=i+1; - beam_matice_id=i+1; //refers to the corresponding material property card - beam_matpar_id=iomodel->numberofnodes-iomodel->numberofnodes2d+1;//refers to the corresponding matpar property card - beam_numpar_id=1; - beam_node_ids[0]=i+1; - beam_node_ids[1]=(int)iomodel->uppernodes[i]; //grid that lays right on top - beam_h[0]=iomodel->thickness[i]; - beam_h[1]=iomodel->thickness[(int)(iomodel->uppernodes[i]-1)]; - beam_s[0]=iomodel->surface[i]; - beam_s[1]=iomodel->surface[(int)(iomodel->uppernodes[i]-1)]; - beam_b[0]=iomodel->bed[i]; - beam_b[1]=iomodel->bed[(int)(iomodel->uppernodes[i]-1)]; - beam_k[0]=iomodel->drag[i]; - beam_k[1]=iomodel->drag[(int)(iomodel->uppernodes[i]-1)]; + /*Create and add penta element to elements dataset: */ + elements->AddObject(new Beam(i,iomodel)); - beam_onbed=(bool)iomodel->gridonbed[i]; + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i,iomodel,2)); - /*Create element properties: */ - beam_properties=new ElementProperties(2,beam_h, beam_s, beam_b, beam_k, NULL,NULL, NULL, UNDEF, NULL, NULL, UNDEF, beam_onbed, (bool)UNDEF, UNDEF, UNDEF, UNDEF); - - /*Create Beam using its constructor:*/ - beam= new Beam(beam_id,beam_node_ids, beam_matice_id, beam_matpar_id, beam_numpar_id, beam_properties); - - /*delete properties: */ - delete beam_properties; - - /*Add tria element to elements dataset: */ - elements->AddObject(beam); - - /*Deal with material property card: */ - matice_mid=i+1; //same as the material id from the geom2 elements. - matice_B=iomodel->B[i]; - matice_n=(double)iomodel->n[1]; //n defined on elements not grids, so take the first value everywhere - - /*Create matice ubeam its constructor:*/ - matice= new Matice(matice_mid,matice_B,matice_n); - - /*Add matice element to materials dataset: */ - materials->AddObject(matice); + } } } - - #ifdef _PARALLEL_ - } //if(my_rank==npart[i]) - #endif } //for (i=0;inumberofnodes;i++) @@ -325,5 +80,4 @@ } //if (strcmp(iomodel->meshtype,"2d")==0) - /*Free data: */ xfree((void**)&iomodel->elements); @@ -337,28 +91,7 @@ xfree((void**)&iomodel->B); xfree((void**)&iomodel->n); - - /*Add one constant material property to materials: */ - matpar_mid=iomodel->numberofnodes-iomodel->numberofnodes2d+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); - - /*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: */ @@ -377,74 +110,18 @@ IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); - - /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); for (i=0;inumberofnodes;i++){ - #ifdef _PARALLEL_ - /*keep only this partition's nodes:*/ - if(my_grids[i]){ - #endif - node_id=i+1; //matlab indexing - + /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ + if(iomodel->my_vertices[i]){ - #ifdef _PARALLEL_ - if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border - node_partitionborder=1; + /*Add vertex to vertices dataset: */ + vertices->AddObject(new Vertex(i,iomodel)); + + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(i,iomodel)); + } - else{ - node_partitionborder=0; - } - #else - node_partitionborder=0; - #endif - - node_x[0]=iomodel->x[i]; - node_x[1]=iomodel->y[i]; - node_x[2]=iomodel->z[i]; - node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); - - node_onbed=(int)iomodel->gridonbed[i]; - node_onsurface=(int)iomodel->gridonsurface[i]; - node_onshelf=(int)iomodel->gridoniceshelf[i]; - node_onsheet=(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,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); - - /*set single point constraints.: */ - if (!iomodel->gridonhutter[i]){ - for(k=1;k<=node_numdofs;k++){ - node->FreezeDof(k); - } - } - - /*Add node to nodes dataset: */ - nodes->AddObject(node); - - #ifdef _PARALLEL_ - } //if(my_grids[i]) - #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(); - materials->Presort(); /*Clean fetched data: */ @@ -461,18 +138,11 @@ xfree((void**)&iomodel->gridonicesheet); xfree((void**)&iomodel->gridoniceshelf); - - /*Keep partitioning information into iomodel*/ - iomodel->epart=epart; - iomodel->my_grids=my_grids; - iomodel->my_bordergrids=my_bordergrids; - - /*Keep partitioning information into iomodel*/ - #ifdef _PARALLEL_ - xfree((void**)&all_numgrids); - xfree((void**)&npart); - VecFree(&gridborder); - #endif - iomodel->npart=npart; + /*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: @@ -481,4 +151,5 @@ *pelements=elements; *pnodes=nodes; + *pvertices=vertices; *pmaterials=materials; Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3422) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3423) @@ -22,11 +22,4 @@ DataSet* vertices = NULL; DataSet* materials = NULL; - - /*Objects: */ - Node* node = NULL; - Vertex* vertex = NULL; - Penta* penta = NULL; - Matice* matice = NULL; - Matpar* matpar = NULL; /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3422) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3423) @@ -98,4 +98,5 @@ /*set single point constraints.: */ + /*FROM DIAGNOSTICHORIZ*/ if (strcmp(iomodel->meshtype,"3d")==0){ /*We have a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ @@ -111,4 +112,5 @@ } } + /*FROM DIAGNOSTICSTOKES*/ /*On a 3d mesh, in stokes formualtions, only stokes grids are free, the others are frozen: */ if (iomodel->borderstokes[i]){ @@ -121,4 +123,10 @@ for(k=1;k<=numdofs;k++){ this->FreezeDof(k); + } + } + /*FROM DIAGNOSTICHUTTER*/ + if (!iomodel->gridonhutter[i]){ + for(k=1;k<=numdofs;k++){ + node->FreezeDof(k); } }