Index: /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp =================================================================== --- /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 3377) +++ /issm/trunk/src/c/ConfigureObjectsx/ConfigureObjectsx.cpp (revision 3378) @@ -18,14 +18,11 @@ extern int my_rank; - _printf_(" Configuring elements...\n"); + //_printf_(" Configuring elements...\n"); elements->Configure(elements,loads,nodes,materials,parameters); - _printf_(" Configuring loads...\n"); - //if (my_rank==2){ - // nodes->Echo(); - //} + //_printf_(" Configuring loads...\n"); loads->Configure(elements,loads,nodes,materials,parameters); - _printf_(" Configuring nodes...\n"); + //_printf_(" Configuring nodes...\n"); nodes->Configure(elements,loads,nodes,materials,parameters); - _printf_(" Configuring parameters...\n"); + //_printf_(" Configuring parameters...\n"); parameters->Configure(elements,loads, nodes, materials,parameters); Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3377) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3378) @@ -2,5 +2,4 @@ * CreateElementsNodesAndMaterialsPrognostic2.c: */ - #include "../../DataSet/DataSet.h" @@ -10,7 +9,7 @@ #include "../../shared/shared.h" #include "../../include/macros.h" +#include "../../include/typedefs.h" #include "../../MeshPartitionx/MeshPartitionx.h" #include "../IoModel.h" - void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ @@ -21,14 +20,14 @@ /*DataSets: */ - DataSet* elements = NULL; - DataSet* nodes = NULL; - DataSet* materials = NULL; + DataSet* elements = NULL; + DataSet* nodes = NULL; + DataSet* materials = NULL; /*Objects: */ - Node* node = NULL; - Tria* tria = NULL; - Penta* penta = NULL; - Matice* matice = NULL; - Matpar* matpar = NULL; + Node* node = NULL; + Tria* tria = NULL; + Penta* penta = NULL; + Matice* matice = NULL; + Matpar* matpar = NULL; /*output: */ @@ -124,4 +123,8 @@ int gridcount; int count; + /*Nodes cloning*/ + double e1,e2; + int i1,i2; + int pos; #endif int first_grid_index; @@ -183,5 +186,4 @@ #endif - /*ids: */ tria_id=i+1; //matlab indexing. @@ -233,8 +235,6 @@ }//for (i=0;ielements); xfree((void**)&iomodel->thickness); xfree((void**)&iomodel->surface); @@ -249,27 +249,58 @@ #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(3*iomodel->numberofelements); - - for (i=0;i<3*iomodel->numberofelements;i++){ - if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + /*If we are in parallel, we must add the nodes that are not in this partition + * but share edges of elements that are in this partition. This is needed for + * the loads as numerical fluxes involve the dofs of all nodes on a given edge*/ + + /*Get edges and elements*/ + IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel_handle,"edges"); + + /*!All elements have been partitioned above, only create elements for this CPU: */ + for (i=0;inumberofedges;i++){ + + /*Get left and right elements*/ + e1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] + e2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] + + /* 1) If the element e1 is in the current partition + * 2) and if the edge of the element is shared by another element + * 3) and if this element is not in the same partition: + * we must clone the nodes on this partition so that the loads (Numericalflux) + * will have access to their properties (dofs,...)*/ + if(my_rank==epart[(int)e1] && !isnan(e2) && my_rank!=epart[(int)e2]){ + + /*1: Get vertices ids*/ + i1=(int)iomodel->edges[4*i+0]; + i2=(int)iomodel->edges[4*i+1]; + + /*2: Get the column where these ids are located in the index*/ + pos==UNDEF; + for(j=0;j<3;j++){ + if (iomodel->elements[3*(int)e2+j]==i1) pos=j+1; + } + ISSMASSERT(pos!=UNDEF); + + /*3: We have the id of the elements and the position of the vertices in the index + * we can now create the corresponding nodes:*/ + my_grids[(int)e2*3+pos-1]=1; + my_grids[(int)e2*3+((pos+1)%3)]=1; } - 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 + } + + /*Free data: */ + xfree((void**)&iomodel->edges); + + /*From the element partitioning, we can determine which nodes are on the inside of this cpu's + *element partition, and which are on its border with other nodes:*/ + gridborder=NewVec(3*iomodel->numberofelements); + + for (i=0;i<3*iomodel->numberofelements;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + VecToMPISerial(&my_bordergrids,gridborder); + #endif @@ -295,5 +326,5 @@ materials->AddObject(matpar); - /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. + /*Ok, let's summarise. Now, every CPU has the following array: my_grids We can therefore determine which grids are internal to this node's partition and which ones are shared with other nodes because they are on the border of this node's partition. Knowing @@ -307,5 +338,4 @@ IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); } - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); @@ -318,5 +348,4 @@ IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); - /*Get number of dofs per node: */ DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); @@ -328,41 +357,49 @@ #ifdef _PARALLEL_ /*!All elements have been partitioned above, only create elements for this CPU: */ - if(my_rank==epart[i]){ + if(my_grids[3*i+j]){ #endif - //Get id of the vertex on which the current node is located - k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing) - ISSMASSERT(k>=0 && knumberofnodes); - - //Get node properties - node_id=i*3+j+1; - node_partitionborder=0; - node_x[0]=iomodel->x[k]; - node_x[1]=iomodel->y[k]; - node_x[2]=iomodel->z[k]; - node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]); - node_onbed=(int)iomodel->gridonbed[k]; - node_onsurface=(int)iomodel->gridonsurface[k]; - node_onshelf=(int)iomodel->gridoniceshelf[k]; - node_onsheet=(int)iomodel->gridonicesheet[k]; - if (strcmp(iomodel->meshtype,"3d")==0){ - if (isnan(iomodel->uppernodes[k])){ - node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + //Get id of the vertex on which the current node is located + k=(int)*(iomodel->elements+elements_width*i+j)-1; //(Matlab to C indexing) + ISSMASSERT(k>=0 && knumberofnodes); + + //Get node properties + node_id=i*3+j+1; + #ifdef _PARALLEL_ + if(my_bordergrids[node_id-1]>1.0) { //this grid belongs to a partition border + node_partitionborder=1; } else{ - node_upper_node_id=(int)iomodel->uppernodes[k]; + node_partitionborder=0; } - } - 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); - - /*Add node to nodes dataset: */ - nodes->AddObject(node); - + #else + node_partitionborder=0; + #endif + node_x[0]=iomodel->x[k]; + node_x[1]=iomodel->y[k]; + node_x[2]=iomodel->z[k]; + node_sigma=(iomodel->z[k]-iomodel->bed[k])/(iomodel->thickness[k]); + node_onbed=(int)iomodel->gridonbed[k]; + node_onsurface=(int)iomodel->gridonsurface[k]; + node_onshelf=(int)iomodel->gridoniceshelf[k]; + node_onsheet=(int)iomodel->gridonicesheet[k]; + if (strcmp(iomodel->meshtype,"3d")==0){ + if (isnan(iomodel->uppernodes[k])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)iomodel->uppernodes[k]; + } + } + 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); + + /*Add node to nodes dataset: */ + nodes->AddObject(node); #ifdef _PARALLEL_ } Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp (revision 3377) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp (revision 3378) @@ -28,5 +28,5 @@ int numericalflux_id; int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES]; - int numericalflux_elem_ids[MAX_NUMERICALFLUX_ELEMS]; + int numericalflux_elem_id; /*Create loads: */ @@ -58,6 +58,5 @@ strcpy(numericalflux_type,"internal"); - numericalflux_elem_ids[0]=(int)e1+1;//id is in matlab index - numericalflux_elem_ids[1]=(int)e2+1;//id is in matlab index + numericalflux_elem_id=(int)e1+1;//id is in matlab index /*Now, we must get the nodes of the 4 nodes located on the edge*/ @@ -85,5 +84,5 @@ strcpy(numericalflux_type,"boundary"); - numericalflux_elem_ids[0]=(int)e1+1; + numericalflux_elem_id=(int)e1+1; /*1: Get vertices ids*/ @@ -104,5 +103,5 @@ } - numericalflux = new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_ids); + numericalflux = new Numericalflux(numericalflux_id,numericalflux_type,numericalflux_node_ids,numericalflux_elem_id); loads->AddObject(numericalflux); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3377) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3378) @@ -220,4 +220,5 @@ idxm=(id-1); value=(double)doflist1[0]; + ISSMASSERT(value>=0); VecSetValues(partition,1,&idxm,&value,INSERT_VALUES); Index: /issm/trunk/src/c/objects/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3377) +++ /issm/trunk/src/c/objects/Numericalflux.cpp (revision 3378) @@ -17,4 +17,6 @@ #include "./objects.h" +extern int my_rank; + /*Object constructors and destructor*/ /*FUNCTION Numericalflux::Numericalflux(){{{1*/ @@ -24,5 +26,5 @@ /*}}}*/ /*FUNCTION Numericalflux::Numericalflux(char numericalflux_type[NUMERICALFLUXSTRING],int numericalflux_fill...){{{1*/ -Numericalflux::Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_ids[MAX_NUMERICALFLUX_ELEMS]){ +Numericalflux::Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_id){ int i; @@ -30,10 +32,8 @@ strcpy(type,numericalflux_type); id=numericalflux_id; - - for(i=0;iEcho(); if (strcmp(type,"internal")==0){ - printf(" element_ids=[%i,%i]\n",element_ids[0],element_ids[1]); - printf(" element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]); - for(i=0;i<2;i++){ - if(elements[i])elements[i]->Echo(); - } printf(" node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]); printf(" node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]); @@ -189,9 +186,4 @@ } else{ - printf(" element_ids=[%i,%i]\n",element_ids[0],element_ids[1]); - printf(" element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]); - for(i=0;i<1;i++){ - if(elements[i])elements[i]->Echo(); - } printf(" node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]); printf(" node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]); @@ -215,13 +207,12 @@ printf(" id: %i\n",id); + printf(" element_id=%i\n",element_id); + printf(" element_offset=%i]\n",element_offset); + if (strcmp(type,"internal")==0){ - printf(" element_ids=[%i,%i]\n",element_ids[0],element_ids[1]); - printf(" element_offsets=[%i,%i]\n",element_offsets[0],element_offsets[1]); printf(" node_ids=[%i,%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2],node_ids[3]); printf(" node_offsets=[%i,%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2],node_offsets[3]); } else{ - printf(" element_ids=%i\n",element_ids[0]); - printf(" element_offsets=%i\n",element_offsets[0]); printf(" node_ids=[%i,%i]\n",node_ids[0],node_ids[1]); printf(" node_offsets=[%i,%i]\n",node_offsets[0],node_offsets[1]); Index: /issm/trunk/src/c/objects/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Numericalflux.h (revision 3377) +++ /issm/trunk/src/c/objects/Numericalflux.h (revision 3378) @@ -24,7 +24,7 @@ /*elements: */ - Element* elements[MAX_NUMERICALFLUX_ELEMS]; - int element_ids[MAX_NUMERICALFLUX_ELEMS]; - int element_offsets[MAX_NUMERICALFLUX_ELEMS]; + Element* element; + int element_id; + int element_offset; /*nodes: */ @@ -36,5 +36,5 @@ Numericalflux(); - Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_ids[MAX_NUMERICALFLUX_ELEMS]); + Numericalflux(int numericalflux_id,char numericalflux_type[NUMERICALFLUXSTRING], int numericalflux_node_ids[MAX_NUMERICALFLUX_NODES],int numericalflux_element_id); Object* copy(); ~Numericalflux(); Index: /issm/trunk/src/c/parallel/prognostic2.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic2.cpp (revision 3377) +++ /issm/trunk/src/c/parallel/prognostic2.cpp (revision 3378) @@ -61,5 +61,5 @@ MPI_Comm_size(MPI_COMM_WORLD,&num_procs); - _printf_("recover , input file name and output file name:\n"); + _printf_("recover input file name and output file name:\n"); inputfilename=argv[2]; outputfilename=argv[3]; Index: /issm/trunk/src/c/parallel/prognostic2_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/prognostic2_core.cpp (revision 3377) +++ /issm/trunk/src/c/parallel/prognostic2_core.cpp (revision 3378) @@ -28,9 +28,8 @@ int numberofdofspernode; int numberofnodes; - int dofs[1]={0}; + int dofs[1]={1}; /*fem prognostic model: */ FemModel* fem_p=NULL; - /*recover fem model: */ Index: /issm/trunk/src/m/classes/public/mesh/meshconvert.m =================================================================== --- /issm/trunk/src/m/classes/public/mesh/meshconvert.m (revision 3377) +++ /issm/trunk/src/m/classes/public/mesh/meshconvert.m (revision 3378) @@ -31,4 +31,5 @@ md.y=bamgmesh_out.Vertices(:,2); md.elements=bamgmesh_out.Triangles(:,1:3); +md.edges=bamgmesh_out.ElementEdges; md.segments=bamgmesh_out.Segments(:,1:3); md.segmentmarkers=bamgmesh_out.Segments(:,4); Index: /issm/trunk/src/m/solutions/jpl/prognostic2_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 3377) +++ /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 3378) @@ -23,5 +23,5 @@ displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']); - results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness'); + %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness'); end %end function