Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 15429) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 15430) @@ -243,4 +243,5 @@ ./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp\ ./modules/ModelProcessorx/NodesPartitioning.cpp\ + ./modules/ModelProcessorx/ProcessMesh2d.cpp\ ./modules/ModelProcessorx/SortDataSets.cpp\ ./modules/ModelProcessorx/UpdateCounters.cpp\ Index: /issm/trunk-jpl/src/c/classes/IoModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 15429) +++ /issm/trunk-jpl/src/c/classes/IoModel.cpp (revision 15430) @@ -64,15 +64,15 @@ /*Initialize permanent data: */ - this->my_elements=NULL; - this->my_nodes=NULL; - this->my_vertices=NULL; - + this->my_elements = NULL; + this->my_nodes = NULL; + this->my_vertices = NULL; FetchData(&this->dim,MeshDimensionEnum); FetchData(&this->numberofvertices,MeshNumberofverticesEnum); FetchData(&this->numberofelements,MeshNumberofelementsEnum); FetchData(&this->elements,NULL,NULL,MeshElementsEnum); - this->edges=NULL; - this->singlenodetoelementconnectivity=NULL; - this->numbernodetoelementconnectivity=NULL; + this->edges = NULL; + this->elementedgeconnectivity = NULL; + this->singlenodetoelementconnectivity = NULL; + this->numbernodetoelementconnectivity = NULL; this->nodecounter=0; @@ -83,6 +83,4 @@ /*FUNCTION IoModel::~IoModel(){{{*/ IoModel::~IoModel(){ - - if(this->constants) delete this->constants; /*Some checks in debugging mode*/ @@ -97,12 +95,17 @@ #endif + if(this->constants) delete this->constants; + xDelete(this->data); xDelete(this->independents); if(this->independent_objects)delete this->independent_objects; + xDelete(this->my_elements); xDelete(this->my_nodes); xDelete(this->my_vertices); + xDelete(this->elements); xDelete(this->edges); + xDelete(this->elementedgeconnectivity); xDelete(this->singlenodetoelementconnectivity); xDelete(this->numbernodetoelementconnectivity); Index: /issm/trunk-jpl/src/c/classes/IoModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/IoModel.h (revision 15429) +++ /issm/trunk-jpl/src/c/classes/IoModel.h (revision 15430) @@ -35,6 +35,8 @@ int numberofvertices; int numberofelements; + int numberofedges; int *elements; int *edges; + int *elementedgeconnectivity; int *singlenodetoelementconnectivity; int *numbernodetoelementconnectivity; Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 15429) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 15430) @@ -33,9 +33,7 @@ /*}}}*/ /*FUNCTION Numericalflux::Numericalflux(int id, int i, IoModel* iomodel, int analysis_type) {{{*/ -Numericalflux::Numericalflux(int numericalflux_id,int i, IoModel* iomodel, int in_analysis_type){ +Numericalflux::Numericalflux(int numericalflux_id,int i,int index,IoModel* iomodel, int in_analysis_type){ /* Intermediary */ - int e1,e2; - int i1,i2; int j; int pos1,pos2,pos3,pos4; @@ -50,19 +48,17 @@ int numericalflux_type; - int numberofelements; - - /*Fetch parameters: */ - iomodel->Constant(&numberofelements,MeshNumberofelementsEnum); - /* Get MatPar id */ - numericalflux_mparid=numberofelements+1; //matlab indexing + numericalflux_mparid=iomodel->numberofelements+1; //matlab indexing + + /*Get edge*/ + int i1 = iomodel->edges[4*index+0]; + int i2 = iomodel->edges[4*index+1]; + int e1 = iomodel->edges[4*index+2]; + int e2 = iomodel->edges[4*index+3]; /*First, see wether this is an internal or boundary edge (if e2=-1)*/ - if (iomodel->Data(MeshEdgesEnum)[4*i+3]==-1.){ //edges are [node1 node2 elem1 elem2] + if(e2==-1){ /* Boundary edge, only one element */ - e1=reCast(iomodel->Data(MeshEdgesEnum)[4*i+2]); - e2=reCast(UNDEF); - num_elems=1; - num_nodes=2; + num_elems=1; num_nodes=2; numericalflux_type=BoundaryEnum; numericalflux_elem_ids[0]=e1; @@ -70,8 +66,5 @@ else{ /* internal edge: connected to 2 elements */ - e1=reCast(iomodel->Data(MeshEdgesEnum)[4*i+2]); - e2=reCast(iomodel->Data(MeshEdgesEnum)[4*i+3]); - num_elems=2; - num_nodes=4; + num_elems=2; num_nodes=4; numericalflux_type=InternalEnum; numericalflux_elem_ids[0]=e1; @@ -80,6 +73,4 @@ /*1: Get vertices ids*/ - i1=reCast(iomodel->Data(MeshEdgesEnum)[4*i+0]); - i2=reCast(iomodel->Data(MeshEdgesEnum)[4*i+1]); numericalflux_vertex_ids[0]=i1; numericalflux_vertex_ids[1]=i2; @@ -93,8 +84,8 @@ pos1=pos2=pos3=pos4=UNDEF; for(j=0;j<3;j++){ - if (iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; - if (iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; - if (iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1; - if (iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1; + if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; + if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; + if(iomodel->elements[3*(e2-1)+j]==i1) pos3=j+1; + if(iomodel->elements[3*(e2-1)+j]==i2) pos4=j+1; } _assert_(pos1!=UNDEF && pos2!=UNDEF && pos3!=UNDEF && pos4!=UNDEF); @@ -112,6 +103,6 @@ pos1=pos2=UNDEF; for(j=0;j<3;j++){ - if (iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; - if (iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; + if(iomodel->elements[3*(e1-1)+j]==i1) pos1=j+1; + if(iomodel->elements[3*(e1-1)+j]==i2) pos2=j+1; } _assert_(pos1!=UNDEF && pos2!=UNDEF); Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 15429) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 15430) @@ -37,5 +37,5 @@ /*Numericalflux constructors,destructors {{{*/ Numericalflux(); - Numericalflux(int numericalflux_id,int i, IoModel* iomodel,int analysis_type); + Numericalflux(int numericalflux_id,int i,int index,IoModel* iomodel,int analysis_type); ~Numericalflux(); /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp (revision 15429) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp (revision 15430) @@ -12,9 +12,7 @@ int element; int stabilization; - int numberofedges; /*Fetch parameters: */ iomodel->Constant(&stabilization,BalancethicknessStabilizationEnum); - iomodel->Constant(&numberofedges,MeshNumberofedgesEnum); /*Recover pointer: */ @@ -28,11 +26,12 @@ /*Get edges and elements*/ - iomodel->FetchData(2,MeshEdgesEnum,ThicknessEnum); + ProcessMesh2d(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel->elements,iomodel->numberofelements,iomodel->numberofvertices); + iomodel->FetchData(1,ThicknessEnum); /*First load data:*/ - for (int i=0;inumberofedges;i++){ /*Get left and right elements*/ - element=reCast(iomodel->Data(MeshEdgesEnum)[4*i+2])-1; //edges are [node1 node2 elem1 elem2] + element=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] /*Now, if this element is not in the partition, pass: */ @@ -40,9 +39,9 @@ /* Add load */ - loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,iomodel,BalancethicknessAnalysisEnum)); + loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,BalancethicknessAnalysisEnum)); } /*Free data: */ - iomodel->DeleteData(2,MeshEdgesEnum,ThicknessEnum); + iomodel->DeleteData(1,ThicknessEnum); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15429) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15430) @@ -49,16 +49,15 @@ if(false){ + /*Quadratic element*/ - int numberofedges; + ProcessMesh2d(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel->elements,iomodel->numberofelements,iomodel->numberofvertices); int element1,element2; bool my_edge; - iomodel->Constant(&numberofedges,MeshNumberofedgesEnum); - - for(int i=0;inumberofedges;i++){ /*Get left and right elements*/ - element1=reCast(iomodel->Data(MeshEdgesEnum)[4*i+2])-1; //edges are [node1 node2 elem1 elem2] - element2=reCast(iomodel->Data(MeshEdgesEnum)[4*i+3])-1; //edges are [node1 node2 elem1 elem2] + element1=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] + element2=iomodel->edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] /*Check whether we should include this edge (element2 is -2 for boundary edges)*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 15429) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 15430) @@ -13,22 +13,22 @@ /*Creation of fem datasets: general drivers*/ -void CreateDataSets(Elements** pelements,Nodes** pnodes,Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,const int nummodels,int analysis_counter); -void CreateElementsVerticesAndMaterials(Elements** pelements,Vertices** pvertices,Materials** pmaterials, IoModel* iomodel,const int nummodels); -void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,int analysis_counter); -void CreateParametersAutodiff(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersDakota(Parameters** pparameters,IoModel* iomodel,char* rootpath,int solution_type,int analysis_type); -void CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersHydrologyDCInefficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); -void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel); +void CreateDataSets(Elements** pelements,Nodes** pnodes,Vertices** pvertices, Materials** pmaterials, Constraints** pconstraints, Loads** ploads,Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,const int nummodels,int analysis_counter); +void CreateElementsVerticesAndMaterials(Elements** pelements,Vertices** pvertices,Materials** pmaterials, IoModel* iomodel,const int nummodels); +void CreateParameters(Parameters** pparameters,IoModel* iomodel,char* rootpath,const int solution_type,int analysis_type,int analysis_counter); +void CreateParametersAutodiff(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); +void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); +void CreateParametersDakota(Parameters** pparameters,IoModel* iomodel,char* rootpath,int solution_type,int analysis_type); +void CreateParametersHydrologyShreve(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); +void CreateParametersHydrologyDCInefficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); +void CreateParametersHydrologyDCEfficient(Parameters** pparameters,IoModel* iomodel,int solution_type,int analysis_type); +void UpdateElementsAndMaterialsControl(Elements* elements,Materials* materials, IoModel* iomodel); /*Creation of fem datasets: specialised drivers: */ /*diagnostic horizontal*/ -void CreateNodesDiagnosticHoriz(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsDiagnosticHoriz(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsDiagnosticHoriz(Loads** ploads, IoModel* iomodel); -void UpdateElementsDiagnosticHoriz(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesDiagnosticHoriz(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsDiagnosticHoriz(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsDiagnosticHoriz(Loads** ploads, IoModel* iomodel); +void UpdateElementsDiagnosticHoriz(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*diagnostic vertical*/ @@ -71,53 +71,54 @@ /*enthalpy:*/ -void CreateNodesEnthalpy(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsEnthalpy(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel); -void UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesEnthalpy(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsEnthalpy(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel); +void UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*hydrology Shreve:*/ -void CreateNodesHydrologyShreve(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsHydrologyShreve(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel); -void UpdateElementsHydrologyShreve(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesHydrologyShreve(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsHydrologyShreve(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsHydrologyShreve(Loads** ploads, IoModel* iomodel); +void UpdateElementsHydrologyShreve(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*hydrology DC:*/ -void CreateNodesHydrologyDCInefficient(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsHydrologyDCInefficient(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsHydrologyDCInefficient(Loads** ploads, IoModel* iomodel); -void UpdateElementsHydrologyDCInefficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); -void CreateNodesHydrologyDCEfficient(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsHydrologyDCEfficient(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsHydrologyDCEfficient(Loads** ploads, IoModel* iomodel); -void UpdateElementsHydrologyDCEfficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesHydrologyDCInefficient(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsHydrologyDCInefficient(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsHydrologyDCInefficient(Loads** ploads, IoModel* iomodel); +void UpdateElementsHydrologyDCInefficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesHydrologyDCEfficient(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsHydrologyDCEfficient(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsHydrologyDCEfficient(Loads** ploads, IoModel* iomodel); +void UpdateElementsHydrologyDCEfficient(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*melting:*/ -void CreateNodesMelting(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsMelting(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsMelting(Loads** ploads, IoModel* iomodel); -void UpdateElementsMelting(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesMelting(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsMelting(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsMelting(Loads** ploads, IoModel* iomodel); +void UpdateElementsMelting(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*prognostic:*/ -void CreateNodesPrognostic(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsPrognostic(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsPrognostic(Loads** ploads, IoModel* iomodel); -void UpdateElementsPrognostic(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesPrognostic(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsPrognostic(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsPrognostic(Loads** ploads, IoModel* iomodel); +void UpdateElementsPrognostic(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*balancedthickness:*/ -void CreateNodesBalancethickness(Nodes** pnodes,IoModel* iomodel); -void CreateConstraintsBalancethickness(Constraints** pconstraints,IoModel* iomodel); -void CreateLoadsBalancethickness(Loads** ploads, IoModel* iomodel); -void UpdateElementsBalancethickness(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); +void CreateNodesBalancethickness(Nodes** pnodes,IoModel* iomodel); +void CreateConstraintsBalancethickness(Constraints** pconstraints,IoModel* iomodel); +void CreateLoadsBalancethickness(Loads** ploads, IoModel* iomodel); +void UpdateElementsBalancethickness(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type); /*transient: */ -void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type); +void UpdateElementsTransient(Elements* elements,Parameters* parameters,IoModel* iomodel,int analysis_counter,int analysis_type); /*partitioning: */ -void ElementsAndVerticesPartitioning(bool** pmy_elements, int** pmy_vertices, IoModel* iomodel); -void NodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices, IoModel* iomodel, bool continuous); +void ElementsAndVerticesPartitioning(bool** pmy_elements, int** pmy_vertices, IoModel* iomodel); +void NodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices, IoModel* iomodel, bool continuous); /*Connectivity*/ -void CreateSingleNodeToElementConnectivity(IoModel* iomodel); -void CreateNumberNodeToElementConnectivity(IoModel* iomodel); +void ProcessMesh2d(int**,int*,int**,int* index,int numberofelements,int numberofvertices); +void CreateSingleNodeToElementConnectivity(IoModel* iomodel); +void CreateNumberNodeToElementConnectivity(IoModel* iomodel); /*Diverse: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp (revision 15429) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/NodesPartitioning.cpp (revision 15430) @@ -20,5 +20,7 @@ void NodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices, IoModel* iomodel, bool continuous){ - /*First thing, this is a new partition for a new analysis_type, therefore, to avoid a leak, erase the nodes partition that might come through pmy_nodes: */ + /*First thing, this is a new partition for a new analysis_type, therefore, + * to avoid a leak, erase the nodes partition that might come through + * pmy_nodes: */ xDelete(*pmy_nodes); @@ -42,18 +44,19 @@ void DiscontinuousGalerkinNodesPartitioning(bool** pmy_nodes,bool* my_elements, int* my_vertices, IoModel* iomodel){ - /*each element has it own nodes (as many as vertices) + additional nodes from neighbouring elements for each edge. This yields to a very different partition for - * the nodes and the vertices. The vertices are similar to continuous galerkin, but the nodes partitioning involves edges, which mess up sorting of - * ids. */ + /* Each element has it own nodes (as many as vertices) + additional nodes + * from neighboring elements for each edge. This yields to a very different + * partition for the nodes and the vertices. The vertices are similar to + * continuous galerkin, but the nodes partitioning involves edges, which + * mess up sorting of ids. */ - /*output: */ - bool* my_nodes=NULL; + /*Intermediaries*/ int i,i1,i2; int cols; int e1,e2; int pos; - int numberofedges; - int *edges = NULL; - int *elements = NULL; + + /*Get edges and elements*/ + ProcessMesh2d(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel->elements,iomodel->numberofelements,iomodel->numberofvertices); /*Build discontinuous node partitioning @@ -67,5 +70,5 @@ /*Allocate*/ - my_nodes=xNewZeroInit(3*iomodel->numberofelements); + bool* my_nodes=xNewZeroInit(3*iomodel->numberofelements); /*First: add all the nodes of all the elements belonging to this cpu*/ @@ -86,13 +89,12 @@ /*Get edges and elements*/ - iomodel->FetchData(&edges,&numberofedges,&cols,MeshEdgesEnum); - if (cols!=4) _error_("field edges should have 4 columns"); + ProcessMesh2d(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel->elements,iomodel->numberofelements,iomodel->numberofvertices); /*!All elements have been partitioned above, only create elements for this CPU: */ - for (i=0;inumberofedges;i++){ /*Get left and right elements*/ - e1=edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] - e2=edges[4*i+3]-1; //edges are [node1 node2 elem1 elem2] + 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 @@ -104,6 +106,6 @@ /*1: Get vertices ids*/ - i1=edges[4*i+0]; - i2=edges[4*i+1]; + i1=iomodel->edges[4*i+0]; + i2=iomodel->edges[4*i+1]; /*2: Get the column where these ids are located in the index*/ @@ -115,5 +117,5 @@ /*3: We have the id of the elements and the position of the vertices in the index * we can now create the corresponding nodes:*/ - if (pos==0){ + if(pos==0){ my_nodes[e2*3+0]=true; my_nodes[e2*3+2]=true; @@ -123,5 +125,5 @@ my_nodes[e2*3+0]=true; } - else if (pos==2){ + else if(pos==2){ my_nodes[e2*3+2]=true; my_nodes[e2*3+1]=true; @@ -133,8 +135,5 @@ } - /*Free data: */ - xDelete(edges); - - /*Assign output pointers:*/ + /*Free data and assign output pointers */ *pmy_nodes=my_nodes; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/ProcessMesh2d.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/ProcessMesh2d.cpp (revision 15430) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/ProcessMesh2d.cpp (revision 15430) @@ -0,0 +1,116 @@ +/*!\file: ProcessMesh2d.cpp + * \brief: create connectivity tables from index + */ + +#include "../../shared/shared.h" + +void ProcessMesh2d(int** pedges,int* pnumberofedges,int** pelement_edge_connectivity,int* index,int numberofelements,int numberofvertices){ + + /*If mesh already processed, return*/ + if(*pedges) return; + + /*Intermediaries*/ + bool exist; + int i,j,v1,v2,v3; + int maxnbe,nbe; + + /*Checks in debugging mode*/ + _assert_(index); + + /*Maximum number of edges*/ + maxnbe = 3*numberofelements; + + /*Initialize intermediaries*/ + int* edgestemp = xNew(maxnbe*4); /*vertex1 vertex2 element1 element2*/ + bool* exchange = xNewZeroInit(maxnbe); + int* element_edge_connectivity = xNew(numberofelements*3); /*edge1 edge2 edge3*/ + for(i=0;i(numberofvertices); + int* next_edge = xNew(maxnbe); + for(i=0;i=0 & v1(head_minv); + xDelete(next_edge); + + /*Create final edges*/ + int* edges = xNew(nbe*4); /*vertex1 vertex2 element1 element2*/ + for(int i=0;i(edgestemp); + xDelete(exchange); + + /*Assign output pointers*/ + if(pedges) *pedges = edges; + else xDelete(edges); + + if(pnumberofedges) *pnumberofedges = nbe; + + if(pelement_edge_connectivity) *pelement_edge_connectivity = element_edge_connectivity; + else xDelete(element_edge_connectivity); +} Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp (revision 15429) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp (revision 15430) @@ -10,15 +10,12 @@ /*Intermediaries*/ - int i; int element; int penpair_ids[2]; int count=0; int stabilization; - int numberofedges; int numvertex_pairing; /*Fetch parameters: */ iomodel->Constant(&stabilization,PrognosticStabilizationEnum); - iomodel->Constant(&numberofedges,MeshNumberofedgesEnum); /*Recover pointer: */ @@ -32,11 +29,12 @@ /*Get edges and elements*/ - iomodel->FetchData(2,MeshEdgesEnum,ThicknessEnum); + ProcessMesh2d(&iomodel->edges,&iomodel->numberofedges,NULL,iomodel->elements,iomodel->numberofelements,iomodel->numberofvertices); + iomodel->FetchData(1,ThicknessEnum); /*First load data:*/ - for (i=0;inumberofedges;i++){ /*Get left and right elements*/ - element=reCast(iomodel->Data(MeshEdgesEnum)[4*i+2])-1; //edges are [node1 node2 elem1 elem2] + element=iomodel->edges[4*i+2]-1; //edges are [node1 node2 elem1 elem2] /*Now, if this element is not in the partition, pass: */ @@ -44,9 +42,9 @@ /* Add load */ - loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,iomodel,PrognosticAnalysisEnum)); + loads->AddObject(new Numericalflux(iomodel->loadcounter+i+1,i,i,iomodel,PrognosticAnalysisEnum)); } /*Free data: */ - iomodel->DeleteData(2,MeshEdgesEnum,ThicknessEnum); + iomodel->DeleteData(1,ThicknessEnum); } @@ -57,5 +55,5 @@ iomodel->FetchData(&nodeonbed,NULL,NULL,MeshVertexonbedEnum); - for(i=0;imy_vertices[reCast(vertex_pairing[2*i+0])-1]){