Index: /issm/trunk/etc/environment_variables.sh =================================================================== --- /issm/trunk/etc/environment_variables.sh (revision 3983) +++ /issm/trunk/etc/environment_variables.sh (revision 3984) @@ -7,5 +7,5 @@ #MATLAB MATLAB_DIR="$ISSM_DIR/externalpackages/matlab/install" -MATLAB_VERSION= +MATLAB_VERSION=R2008a #MPI Index: /issm/trunk/externalpackages/matlab/install.sh =================================================================== --- /issm/trunk/externalpackages/matlab/install.sh (revision 3983) +++ /issm/trunk/externalpackages/matlab/install.sh (revision 3984) @@ -8,13 +8,13 @@ # MAC symlink matlab to root matlab -#ln -s /Applications/MATLAB_R2008a/ install +ln -s /Applications/MATLAB_R2008a/ install #export MATLAB_VERSION=R2008a # LINUX symlink matlab to root matlab -ln -s /usr/local/pkgs/matlab-7.6/ install +#ln -s /usr/local/pkgs/matlab-7.6/ install #initialize MATLAB_VERSION in etc/environment_variables -cat $ISSM_DIR/etc/environment_variables.sh | sed -e "/MATLAB_VERSION/d" | sed "/^MATLAB_DIR.*$/a\MATLAB_VERSION=$MATLAB_VERSION" > $ISSM_DIR/etc/environment_variables.sh.bak -mv $ISSM_DIR/etc/environment_variables.sh.bak $ISSM_DIR/etc/environment_variables.sh -cat $ISSM_DIR/etc/environment_variables.csh | sed -e "/MATLAB_VERSION/d" | sed "/^set MATLAB_DIR.*$/a\set MATLAB_VERSION=$MATLAB_VERSION" > $ISSM_DIR/etc/environment_variables.csh.bak -mv $ISSM_DIR/etc/environment_variables.csh.bak $ISSM_DIR/etc/environment_variables.csh +#cat $ISSM_DIR/etc/environment_variables.sh | sed -e "/MATLAB_VERSION/d" | sed "/^MATLAB_DIR.*$/a\MATLAB_VERSION=$MATLAB_VERSION" > $ISSM_DIR/etc/environment_variables.sh.bak +#mv $ISSM_DIR/etc/environment_variables.sh.bak $ISSM_DIR/etc/environment_variables.sh +#cat $ISSM_DIR/etc/environment_variables.csh | sed -e "/MATLAB_VERSION/d" | sed "/^set MATLAB_DIR.*$/a\set MATLAB_VERSION=$MATLAB_VERSION" > $ISSM_DIR/etc/environment_variables.csh.bak +#mv $ISSM_DIR/etc/environment_variables.csh.bak $ISSM_DIR/etc/environment_variables.csh Index: /issm/trunk/externalpackages/metis/configs/macosx64/Makefile.in.patch =================================================================== --- /issm/trunk/externalpackages/metis/configs/macosx64/Makefile.in.patch (revision 3983) +++ /issm/trunk/externalpackages/metis/configs/macosx64/Makefile.in.patch (revision 3984) @@ -8,3 +8,3 @@ < OPTFLAGS = -DLINUX --- -> OPTFLAGS = -O2 +> OPTFLAGS = -O2 -m64 Index: /issm/trunk/externalpackages/triangle/configs/macosx64/configure.make =================================================================== --- /issm/trunk/externalpackages/triangle/configs/macosx64/configure.make (revision 3983) +++ /issm/trunk/externalpackages/triangle/configs/macosx64/configure.make (revision 3984) @@ -1,3 +1,3 @@ -CC=gcc +CC=gcc -m64 CSWITCHES = $(CFLAGS) -I/usr/X11R6/include -L/usr/X11R6/lib -I$(MATLAB_DIR)/extern/include -fPIC -I$(MATLAB_DIR)/include TRILIBDEFS = -DTRILIBRARY Index: /issm/trunk/m4/matlab.m4 =================================================================== --- /issm/trunk/m4/matlab.m4 (revision 3983) +++ /issm/trunk/m4/matlab.m4 (revision 3984) @@ -312,4 +312,10 @@ MATLAB_MINOR=`echo $MATLAB_VERSION | sed -e 's/^@<:@0-9@:>@*\.\(@<:@0-9@:>@*\).*/\1/'` ;; + R2010a) + MATLAB_VERSION=R2010a + MATLAB_MAJOR=7 + MATLAB_MINOR=8 + ;; + R2009a) MATLAB_VERSION=R2009a Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3983) +++ /issm/trunk/src/c/Makefile.am (revision 3984) @@ -298,4 +298,6 @@ ./EnumDefinitions/StringAsEnum.cpp\ ./modules/ModelProcessorx/ModelProcessorx.h\ + ./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp + ./modules/ModelProcessorx/NodesPartitioning.cpp ./modules/ModelProcessorx/Partitioning.cpp\ ./modules/ModelProcessorx/CreateDataSets.cpp\ @@ -748,4 +750,6 @@ ./EnumDefinitions/StringAsEnum.cpp\ ./modules/ModelProcessorx/ModelProcessorx.h\ + ./modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp + ./modules/ModelProcessorx/NodesPartitioning.cpp ./modules/ModelProcessorx/Partitioning.cpp\ ./modules/ModelProcessorx/CreateDataSets.cpp\ Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp (revision 3984) @@ -34,5 +34,5 @@ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ /*Fetch data needed: */ @@ -75,5 +75,5 @@ } - else{ // if (strcmp(meshtype,"2d")==0) + else{ // if (dim==2) /*Fetch data needed: */ @@ -122,5 +122,5 @@ - } //if (strcmp(meshtype,"2d")==0) + } //if (dim==2) /*Add new constrant material property to materials, at the end: */ @@ -128,5 +128,5 @@ /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedthickness2/CreateElementsNodesAndMaterialsBalancedthickness2.cpp (revision 3984) @@ -37,5 +37,5 @@ /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ /*Fetch data needed: */ @@ -77,7 +77,7 @@ } - else{ // if (strcmp(meshtype,"2d")==0) + else{ // if (dim==2) ISSMERROR("not implemented yet"); - } //if (strcmp(meshtype,"2d")==0) + } //if (dim==2) /*Add new constrant material property tgo materials, at the end: */ @@ -95,5 +95,5 @@ IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (revision 3984) @@ -34,5 +34,5 @@ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ /*Fetch data needed: */ @@ -75,5 +75,5 @@ } - else{ // if (strcmp(meshtype,"2d")==0) + else{ // if (dim==2) /*Fetch data needed: */ @@ -122,5 +122,5 @@ - } //if (strcmp(meshtype,"2d")==0) + } //if (dim==2) /*Add new constrant material property to materials, at the end: */ @@ -128,5 +128,5 @@ /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 3984) @@ -16,92 +16,106 @@ -void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_type){ +void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_type,int nummodels,int analysis_counter){ - /*First, partition elements and vertices. Nodes will partitioned on a per analysis_type basis: */ - VerticesPartitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_bordervertices, iomodel, iomodel_handle); + bool continuous=true; + DataSet* elements=NULL; + + /*Create elements, vertices and materials, independent of analysis_type: */ + CreateElementsVerticesAndMaterials(pelements, pvertices, pmaterials, iomodel,iomodel_handle,nummodels); - /*create parameters common to all solutions: */ + /*Recover elements, for future update: */ + elements=*pelements; + + /*Now, branch onto analysis dependent model generation: */ + switch(analysis_type){ + case DiagnosticHorizAnalysisEnum: + CreateNodesDiagnosticHoriz(pnodes, iomodel,iomodel_handle); + CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle); + UpdateElementsDiagnosticHoriz(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case DiagnosticVertAnalysisEnum: + CreateNodesDiagnosticVert(pnodes, iomodel,iomodel_handle); + CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle); + UpdateElementsDiagnosticVert(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case DiagnosticStokesAnalysisEnum: + CreateNodesDiagnosticStokes(pnodes, iomodel,iomodel_handle); + CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle); + UpdateElementsDiagnosticStokes(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case DiagnosticHutterAnalysisEnum: + CreateNodesDiagnosticHutter(pnodes, iomodel,iomodel_handle); + CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle); + UpdateElementsDiagnosticHutter(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case SlopecomputeAnalysisEnum: + CreateNodesSlopeCompute(pnodes, iomodel,iomodel_handle); + CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle); + CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle); + UpdateElementsSlopeCompute(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case ThermalAnalysisEnum: + CreateNodesThermal(pnodes, iomodel,iomodel_handle); + CreateConstraintsThermal(pconstraints,iomodel,iomodel_handle); + CreateLoadsThermal(ploads,iomodel,iomodel_handle); + UpdateElementsThermal(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case MeltingAnalysisEnum: + CreateNodesMelting(pnodes, iomodel,iomodel_handle); + CreateConstraintsMelting(pconstraints,iomodel,iomodel_handle); + CreateLoadsMelting(ploads,iomodel,iomodel_handle); + UpdateElementsMelting(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case PrognosticAnalysisEnum: + CreateNodesPrognostic(pnodes, iomodel,iomodel_handle); + CreateConstraintsPrognostic(pconstraints,iomodel,iomodel_handle); + CreateLoadsPrognostic(ploads,iomodel,iomodel_handle); + UpdateElementsPrognostic(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case Prognostic2AnalysisEnum: + CreateNodesPrognostic2(pnodes, iomodel,iomodel_handle); + CreateConstraintsPrognostic2(pconstraints,iomodel,iomodel_handle); + CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle); + UpdateElementsPrognostic2(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case BalancedthicknessAnalysisEnum: + CreateNodesBalancedthickness(pnodes, iomodel,iomodel_handle); + CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle); + CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle); + UpdateElementsBalancedthickness(elements,iomodel,iomodel_handle,analysis_counter); + break; + + case Balancedthickness2AnalysisEnum: + CreateNodesBalancedthickness2(pnodes, iomodel,iomodel_handle); + CreateConstraintsBalancedthickness2(pconstraints,iomodel,iomodel_handle); + CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle); + UpdateElementsBalancedthickness2(elements,iomodel,iomodel_handle,analysis_counter); + break; + case BalancedvelocitiesAnalysisEnum: + CreateNodesBalancedvelocities(pnodes, iomodel,iomodel_handle); + CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle); + CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle); + UpdateElementsBalancedvelocities(elements,iomodel,iomodel_handle,analysis_counter); + break; + default: + ISSMERROR("%s%s%s"," analysis_type: ",EnumAsString(analysis_type)," not supported yet!"); + } + + /*Generate objects that are not dependent on any analysis_type: */ CreateParameters(pparameters,iomodel,iomodel_handle); - CreateParametersControl(pparameters,iomodel,iomodel_handle); - - /*This is just a high level driver: */ - if (analysis_type==DiagnosticHorizAnalysisEnum){ - CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==DiagnosticVertAnalysisEnum){ - - CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==DiagnosticStokesAnalysisEnum){ - - CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==DiagnosticHutterAnalysisEnum){ - - CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==SlopecomputeAnalysisEnum){ - - CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes, pvertices,pmaterials, iomodel,iomodel_handle); - CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle); - CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==ThermalAnalysisEnum){ - - CreateElementsNodesAndMaterialsThermal(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsThermal(pconstraints,iomodel,iomodel_handle); - CreateLoadsThermal(ploads,iomodel,iomodel_handle); - - } - else if (analysis_type==MeltingAnalysisEnum){ - - CreateElementsNodesAndMaterialsMelting(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsMelting(pconstraints,iomodel,iomodel_handle); - CreateLoadsMelting(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==PrognosticAnalysisEnum){ - - CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsPrognostic(pconstraints,iomodel,iomodel_handle); - CreateLoadsPrognostic(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==Prognostic2AnalysisEnum){ - - CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsPrognostic2(pconstraints,iomodel,iomodel_handle); - CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==BalancedthicknessAnalysisEnum){ - - CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsBalancedthickness(pconstraints,iomodel,iomodel_handle); - CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==Balancedthickness2AnalysisEnum){ - - CreateElementsNodesAndMaterialsBalancedthickness2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsBalancedthickness2(pconstraints,iomodel,iomodel_handle); - CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle); - } - else if (analysis_type==BalancedvelocitiesAnalysisEnum){ - - CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pvertices,pmaterials, iomodel,iomodel_handle); - CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle); - CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle); - } - else{ - ISSMERROR("%s%i%s%i%s"," analysis_type: ",analysis_type," analysis_type: ",analysis_type," not supported yet!"); - - } CreateParametersQmu(pparameters,iomodel,iomodel_handle); - } Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 3984) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 3984) @@ -0,0 +1,103 @@ +/* + * CreateElementsNodesAndMaterialsDiagnosticHoriz.c: + */ + +#include "../../../DataSet/DataSet.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateElementsVerticesAndMaterials(DataSet** pelements,DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle,int nummodels){ + + /*Intermediary*/ + int i,j,k,n; + + /*DataSets: */ + DataSet* elements = NULL; + DataSet* vertices = NULL; + DataSet* materials = NULL; + + /*Did we already create the elements? : */ + if(*pelements)return; + + /*First create the elements, vertices, nodes and material properties, if they don't already exist */ + elements = new DataSet(ElementsEnum); + vertices = new DataSet(VerticesEnum); + materials = new DataSet(MaterialsEnum); + + /*Partition elements and vertices and nodes: */ + Partitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_nodes, &iomodel->my_bordervertices, iomodel, iomodel_handle); + + /*First, partition elements and vertices. Nodes will partitioned on a per analysis_type basis. If partitining already done, ignore: */ + ElementsAndVerticesPartitioning(&iomodel->my_elements, &iomodel->my_vertices, &iomodel->my_bordervertices, iomodel, iomodel_handle); + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements"); + IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements"); + + /*Create elements and materials: */ + for (i=0;inumberofelements;i++){ + + if(iomodel->my_elements[i]){ + + /*Create and add tria element to elements dataset: */ + if(iomodel->dim==2)elements->AddObject(new Tria(i+1,i,iomodel,nummodels)); + else elements->AddObject(new Penta(i+1,i,iomodel,nummodels)); + + /*Create and add material property to materials dataset: */ + materials->AddObject(new Matice(i+1,i,iomodel)); + } + } + + /*Free data: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->upperelements); + xfree((void**)&iomodel->lowerelements); + + /*Add new constrant material property tgo materials, at the end: */ + materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials + + /*Create vertices: */ + IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); + IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); + IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + + for (i=0;inumberofvertices;i++){ + + /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ + if(iomodel->my_vertices[i]){ + + /*Add vertex to vertices dataset: */ + vertices->AddObject(new Vertex(i+1,i,iomodel)); + } + } + + /*Clean fetched data: */ + xfree((void**)&iomodel->x); + xfree((void**)&iomodel->y); + xfree((void**)&iomodel->z); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->thickness); + + /*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(); + vertices->Presort(); + materials->Presort(); + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pvertices=vertices; + *pmaterials=materials; + +} Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateNumberNodeToElementConnectivity.cpp (revision 3984) @@ -33,5 +33,5 @@ /*Get element width (3 or 6)*/ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ elementswidth=3; } Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 3984) @@ -28,5 +28,5 @@ parameters = new Parameters(ParametersEnum); - if (strcmp(iomodel->meshtype,"2d")==0) parameters->AddObject(new IntParam(DimEnum,2)); + if (iomodel->dim==2) parameters->AddObject(new IntParam(DimEnum,2)); else parameters->AddObject(new IntParam(DimEnum,3)); parameters->AddObject(new IntParam(AnalysisTypeEnum,iomodel->analysis_type)); Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateSingleNodeToElementConnectivity.cpp (revision 3984) @@ -34,5 +34,5 @@ /*Get element width (3 or 6)*/ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ elementswidth=3; } Index: sm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 3983) +++ (revision ) @@ -1,234 +1,0 @@ -/* - * CreateElementsNodesAndMaterialsDiagnosticHoriz.c: - */ - -#include "../../../DataSet/DataSet.h" -#include "../../../toolkits/toolkits.h" -#include "../../../io/io.h" -#include "../../../EnumDefinitions/EnumDefinitions.h" -#include "../../../objects/objects.h" -#include "../../../shared/shared.h" -#include "../../MeshPartitionx/MeshPartitionx.h" -#include "../../../include/include.h" -#include "../ModelProcessorx.h" - -void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ - - /*Intermediary*/ - int i,j,k,n; - - /*DataSets: */ - DataSet* elements = NULL; - DataSet* nodes = NULL; - DataSet* vertices = NULL; - DataSet* materials = NULL; - - /*First create the elements, nodes and material properties: */ - elements = new DataSet(ElementsEnum); - nodes = new DataSet(NodesEnum); - vertices = new DataSet(VerticesEnum); - materials = new DataSet(MaterialsEnum); - - /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if (!iomodel->ismacayealpattyn)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); - - /*2d mesh: */ - if (strcmp(iomodel->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_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); - IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); - IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); - IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); - IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); - IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); - IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); - IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); - IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); - IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); - if (iomodel->control_analysis){ - IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs"); - IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); - IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); - } - - for (i=0;inumberofelements;i++){ - - if(iomodel->my_elements[i]){ - - if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ //elements of type 1 are Hutter type Tria. Don't create this elements. - - /*Create and add tria element to elements dataset: */ - elements->AddObject(new Tria(i+1,i,iomodel)); - - /*Create and add material property to materials dataset: */ - materials->AddObject(new Matice(i+1,i,iomodel,3)); - } - } - }//for (i=0;ielements); - xfree((void**)&iomodel->elements_type); - xfree((void**)&iomodel->thickness); - xfree((void**)&iomodel->surface); - xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->drag_coefficient); - xfree((void**)&iomodel->drag_p); - xfree((void**)&iomodel->drag_q); - xfree((void**)&iomodel->rheology_B); - xfree((void**)&iomodel->rheology_n); - xfree((void**)&iomodel->elementoniceshelf); - xfree((void**)&iomodel->elementonwater); - xfree((void**)&iomodel->vx); - xfree((void**)&iomodel->vy); - if (iomodel->control_analysis){ - xfree((void**)&iomodel->vx_obs); - xfree((void**)&iomodel->vy_obs); - xfree((void**)&iomodel->weights); - } - - } - else{ // if (strcmp(meshtype,"2d")==0) - - /*Fetch data needed: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); - 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_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); - IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); - IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); - IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); - IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); - IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); - 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->vx,NULL,NULL,iomodel_handle,"vx"); - IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); - IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz"); - if (iomodel->control_analysis){ - IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs"); - IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); - IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); - } - IoModelFetchData(&iomodel->upperelements,NULL,NULL,iomodel_handle,"upperelements"); - IoModelFetchData(&iomodel->lowerelements,NULL,NULL,iomodel_handle,"lowerelements"); - - - for (i=0;inumberofelements;i++){ - if(iomodel->my_elements[i]){ - if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum || *(iomodel->elements_type+2*i+0)==PattynFormulationEnum){ //elements of type 1 are Hutter type Tria. Don't create this elements. - /*Create and add penta element to elements dataset: */ - elements->AddObject(new Penta(i+1,i,iomodel)); - - /*Create and add material property to materials dataset: */ - materials->AddObject(new Matice(i+1,i,iomodel,6)); - - } - }//if(my_elements[i]) - - }//for (i=0;ielements); - xfree((void**)&iomodel->elements_type); - xfree((void**)&iomodel->thickness); - xfree((void**)&iomodel->surface); - xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->drag_coefficient); - xfree((void**)&iomodel->drag_p); - xfree((void**)&iomodel->drag_q); - xfree((void**)&iomodel->rheology_n); - xfree((void**)&iomodel->rheology_B); - xfree((void**)&iomodel->elementoniceshelf); - xfree((void**)&iomodel->elementonbed); - xfree((void**)&iomodel->elementonsurface); - xfree((void**)&iomodel->elementonwater); - xfree((void**)&iomodel->vx); - xfree((void**)&iomodel->vy); - xfree((void**)&iomodel->vz); - if (iomodel->control_analysis){ - xfree((void**)&iomodel->vx_obs); - xfree((void**)&iomodel->vy_obs); - xfree((void**)&iomodel->weights); - } - xfree((void**)&iomodel->upperelements); - xfree((void**)&iomodel->lowerelements); - - - - } //if (strcmp(meshtype,"2d")==0) - - /*Add new constrant material property tgo materials, at the end: */ - materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel));//put it at the end of the materials - - /*Create nodes and vertices: */ - IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); - IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); - IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z"); - IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); - IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); - IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); - IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); - IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); - IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); - IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); - if (strcmp(iomodel->meshtype,"3d")==0){ - IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); - IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); - } - - for (i=0;inumberofvertices;i++){ - - /*vertices and nodes (same number, as we are running continuous galerkin formulation: */ - if(iomodel->my_vertices[i]){ - - /*Add vertex to vertices dataset: */ - vertices->AddObject(new Vertex(i+1,i,iomodel)); - - /*Add node to nodes dataset: */ - nodes->AddObject(new Node(i+1,i,iomodel)); - - } - } - - /*Clean fetched data: */ - xfree((void**)&iomodel->deadgrids); - xfree((void**)&iomodel->x); - xfree((void**)&iomodel->y); - xfree((void**)&iomodel->z); - xfree((void**)&iomodel->thickness); - xfree((void**)&iomodel->bed); - xfree((void**)&iomodel->gridonbed); - xfree((void**)&iomodel->gridonsurface); - xfree((void**)&iomodel->gridonhutter); - xfree((void**)&iomodel->uppernodes); - xfree((void**)&iomodel->gridonicesheet); - xfree((void**)&iomodel->gridoniceshelf); - - /*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: - - /*Assign output pointer: */ - *pelements=elements; - *pnodes=nodes; - *pvertices=vertices; - *pmaterials=materials; - -} Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 3984) @@ -39,5 +39,5 @@ /*Retrieve element to which this icefront belongs: */ - if (strcmp(iomodel->meshtype,"2d")==0) segment_width=4; + if (iomodel->dim==2) segment_width=4; else segment_width=6; element=(int)(*(iomodel->pressureload+segment_width*i+segment_width-2)-1); //element is in the penultimate column (grid1 grid2 ... elem fill) Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 3984) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 3984) @@ -0,0 +1,69 @@ +/* + * CreateNodesDiagnosticHoriz.c: + */ + +#include "../../../DataSet/DataSet.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateNodesDiagnosticHoriz(DataSet** pnodes, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + /*Intermediary*/ + int i; + bool galerkin_continuous=true; + + /*DataSets: */ + DataSet* nodes = NULL; + + /*First create the elements, vertices, nodes and material properties, if they don't already exist */ + nodes = new DataSet(NodesEnum); + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!iomodel->ismacayealpattyn)goto cleanup_and_return; + + /*Galerkin continus partition of nodes: */ + NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements, iomodel->my_vertices, iomodel->my_bordervertices, iomodel, iomodel_handle,galerkin_continuous); + + /*Create nodes: */ + IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); + IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); + IoModelFetchData(&iomodel->gridonhutter,NULL,NULL,iomodel_handle,"gridonhutter"); + IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); + IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + if (iomodel->dim==3){ + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); + } + + for (i=0;inumberofvertices;i++){ + + if(iomodel->my_vertices[i]){ + + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(i+1,i,iomodel,DiagnosticHorizAnalysisEnum)); + } + } + + /*Clean fetched data: */ + xfree((void**)&iomodel->gridonbed); + xfree((void**)&iomodel->gridonsurface); + xfree((void**)&iomodel->gridonhutter); + xfree((void**)&iomodel->gridonicesheet); + xfree((void**)&iomodel->gridoniceshelf); + xfree((void**)&iomodel->deadgrids); + xfree((void**)&iomodel->uppernodes); + + /*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: */ + nodes->Presort(); + + /*Assign output pointer: */ + *pnodes=nodes; + +} Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 3984) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 3984) @@ -0,0 +1,83 @@ +/* + * UpdateElementsDiagnosticHoriz: + */ + +#include "../../../DataSet/DataSet.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsDiagnosticHoriz(DataSet* elements, IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_counter){ + + /*Intermediary*/ + int i; + Element* element=NULL; + + /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ + if (!iomodel->ismacayealpattyn)goto cleanup_and_return; + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); + 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_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); + IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); + IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); + IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); + IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); + IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); + IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); + if (iomodel->control_analysis){ + IoModelFetchData(&iomodel->vx_obs,NULL,NULL,iomodel_handle,"vx_obs"); + IoModelFetchData(&iomodel->vy_obs,NULL,NULL,iomodel_handle,"vy_obs"); + IoModelFetchData(&iomodel->weights,NULL,NULL,iomodel_handle,"weights"); + } + + if (iomodel->dim==3){ + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz"); + } + + /*Update elements: */ + for (i=0;iSize();i++){ + + element=(Element*)elements->GetObjectByOffset(i); + element->Update(iomodel,analysis_counter); + + } + + /*Free data: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements_type); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->drag_p); + xfree((void**)&iomodel->drag_q); + xfree((void**)&iomodel->rheology_n); + xfree((void**)&iomodel->rheology_B); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elementonwater); + xfree((void**)&iomodel->vx); + xfree((void**)&iomodel->vy); + xfree((void**)&iomodel->vz); + if (iomodel->control_analysis){ + xfree((void**)&iomodel->vx_obs); + xfree((void**)&iomodel->vy_obs); + xfree((void**)&iomodel->weights); + } + xfree((void**)&iomodel->upperelements); + xfree((void**)&iomodel->lowerelements); +} Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticHutter/CreateElementsNodesAndMaterialsDiagnosticHutter.cpp (revision 3984) @@ -59,5 +59,5 @@ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ for (i=0;inumberofelements;i++){ @@ -75,5 +75,5 @@ } } //for (i=0;inumberofvertices;i++) - } //if (strcmp(iomodel->meshtype,"2d")==0) + } //if (iomodel->dim==2) else{ @@ -96,5 +96,5 @@ } //for (i=0;inumberofelements;i++) - } //if (strcmp(iomodel->meshtype,"2d")==0) + } //if (iomodel->dim==2) /*Free data: */ @@ -117,5 +117,5 @@ /*Add new constrant material property to materials, at the end: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ materials->AddObject(new Matpar(iomodel->numberofelements+1,iomodel)); //put it at the end of the materials } @@ -125,5 +125,5 @@ /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticStokes/CreateElementsNodesAndMaterialsDiagnosticStokes.cpp (revision 3984) @@ -31,5 +31,5 @@ /*Now, do we have Stokes elements?*/ - if (strcmp(iomodel->meshtype,"2d")==0) goto cleanup_and_return; + if (iomodel->dim==2) goto cleanup_and_return; if (!iomodel->isstokes) goto cleanup_and_return; Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateConstraintsDiagnosticVert.cpp (revision 3984) @@ -21,5 +21,5 @@ /*return if 2d mesh*/ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; + if (iomodel->dim==2)goto cleanup_and_return; /*Fetch data: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/DiagnosticVert/CreateElementsNodesAndMaterialsDiagnosticVert.cpp (revision 3984) @@ -31,5 +31,5 @@ /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; + if (iomodel->dim==2)goto cleanup_and_return; /*Partition elements and vertices and nodes: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (revision 3984) +++ /issm/trunk/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp (revision 3984) @@ -0,0 +1,147 @@ +/*!\file: ElementsAndVerticesPartitioning.cpp + * \brief: partition elements and nodes and vertices + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../io/io.h" +#include "../../include/include.h" +#include "../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx/ModelProcessorx.h" + +void ElementsAndVerticesPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ + + int i; + + extern int my_rank; + extern int num_procs; + + /*output: */ + bool* my_elements=NULL; + bool* my_vertices=NULL; + bool* my_bordervertices=NULL; + + /*intermediary: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int elements_width; //number of columns in elements (2d->3, 3d->6) + Vec bordervertices=NULL; + double* serial_bordervertices=NULL; + int el1,el2; + + /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */ + if (*pmy_elements)return; + + /*Number of vertices per elements, needed to correctly retrieve data: */ + if(iomodel->dim==2) elements_width=3; //tria elements + else elements_width=6; //penta elements + + #ifdef _PARALLEL_ + /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ + if(iomodel->dim==2){ + /*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->numberofvertices,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofvertices2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->dim,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements2d); + + #else + /*In serial mode, epart is full of 0: all elements belong to cpu 0: */ + epart=(int*)xcalloc(iomodel->numberofelements,sizeof(int)); + #endif + + /*Deal with rifts, they have to be included into one partition only, not several: */ + IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo"); + + for(i=0;inumrifts;i++){ + el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing + el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing + epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; + } + + /*Free rifts: */ + xfree((void**)&iomodel->riftinfo); + + /*Used later on: */ + my_vertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); + my_elements=(bool*)xcalloc(iomodel->numberofelements,sizeof(bool)); + + /*Start figuring out, out of the partition, which elements belong to this cpu: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + for (i=0;inumberofelements;i++){ + + /*!All elements have been partitioned above, only deal with elements for this cpu: */ + if(my_rank==epart[i]){ + + my_elements[i]=1; + + /*Now that we are here, we can also start building the list of vertices belonging to this cpu 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 vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices + will hold which vertices belong to this partition*/ + my_vertices[(int)*(iomodel->elements+elements_width*i+0)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+1)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + + if(elements_width==6){ + my_vertices[(int)*(iomodel->elements+elements_width*i+3)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+4)-1]=1; + my_vertices[(int)*(iomodel->elements+elements_width*i+5)-1]=1; + } + } + }//for (i=0;ielements); + + #ifdef _PARALLEL_ + /*From the element partitioning, we can determine which vertices are on the inside of this cpu's + *element partition, and which are on its border with other vertices:*/ + bordervertices=NewVec(iomodel->numberofvertices); + + for (i=0;inumberofvertices;i++){ + if(my_vertices[i])VecSetValue(bordervertices,i,1,ADD_VALUES); + } + VecAssemblyBegin(bordervertices); + VecAssemblyEnd(bordervertices); + + VecToMPISerial(&serial_bordervertices,bordervertices); + + /*now go through serial_bordervertices, and booleanize it: */ + my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); + for(i=0;inumberofvertices;i++){ + if(serial_bordervertices[i]>1)my_bordervertices[i]=1; + } + #else + /*No border vertices: */ + my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); + #endif + + /*Free ressources:*/ + xfree((void**)&npart); + xfree((void**)&epart); + xfree((void**)&serial_bordervertices); + VecFree(&bordervertices); + + /*Assign output pointers:*/ + *pmy_elements=my_elements; + *pmy_vertices=my_vertices; + *pmy_bordervertices=my_bordervertices; +} + + Index: /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateElementsNodesAndMaterialsMelting.cpp (revision 3984) @@ -92,5 +92,5 @@ /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Melting/CreateLoadsMelting.cpp (revision 3984) @@ -18,5 +18,5 @@ /*if 2d: return*/ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; + if (iomodel->dim==2)goto cleanup_and_return; /*Create loads: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 3984) @@ -13,84 +13,91 @@ #include "../../io/io.h" +/*Creation of fem datasets: general drivers*/ +void CreateDataSets(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_type,int nummodels,int analysis_counter); +void CreateElementsVerticesAndMaterials(DataSet** pelements,DataSet** pvertices,DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle,int nummodels); +void CreateParameters(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); -/*Creation of fem datasets: general drivers*/ -void CreateDataSets(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_type); -void CreateParameters(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); - -/*Create of fem datasets: specialised drivers: */ +/*Creation of fem datasets: specialised drivers: */ /*diagnostic horizontal*/ -void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesDiagnosticHoriz(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElementsDiagnosticHoriz(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*diagnostic vertical*/ -void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesDiagnosticVert(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*diagnostic hutter*/ -void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesDiagnosticHutter(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*diagnostic stokes*/ -void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesDiagnosticStokes(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*slope compute*/ -void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesSlopeCompute(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - -/*control:*/ -void CreateParametersControl(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*thermal:*/ -void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesThermal(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*melting:*/ -void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesMelting(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*prognostic:*/ -void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesPrognostic(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*prognostic2:*/ -void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesPrognostic2(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsPrognostic2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*balancedthickness:*/ -void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesBalancedthickness(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); -void CreateElementsNodesAndMaterialsBalancedthickness2(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesBalancedthickness2(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsBalancedthickness2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - +void CreateLoadsBalancedthickness2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*balancedvelocities:*/ -void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes,DataSet** pvertices, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateNodesBalancedvelocities(pnodes,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); void CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); -void CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - -/*qmu: */ -void CreateParametersQmu(Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); +void CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); +void UpdateElements(DataSet* elements,IoModel* iomodel_handle,ConstDataHandle iomodel_handle); /*partitioning: */ -void VerticesPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); -void NodesPartitioning(bool** pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle,bool continuous); +void ElementsAndVerticesPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); +void NodesPartitioning(bool** pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle,bool continuous); -/*Conectivity*/ -void CreateSingleNodeToElementConnectivity(IoModel* iomodel); -void CreateNumberNodeToElementConnectivity(IoModel* iomodel); +/*Connectivity*/ +void CreateSingleNodeToElementConnectivity(IoModel* iomodel); +void CreateNumberNodeToElementConnectivity(IoModel* iomodel); #endif Index: /issm/trunk/src/c/modules/ModelProcessorx/NodesPartitioning.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/NodesPartitioning.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/NodesPartitioning.cpp (revision 3984) @@ -23,4 +23,8 @@ void NodesPartitioning(bool** pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle,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: */ + xfree((void**)pmy_nodes); + + /*Now, depending on whether we are running galerkin discontinous or continuous elements, carry out a different partition of the nodes: */ if(continuous==true){ ContinuousGalerkinNodesPartitioning(pmy_nodes,my_elements, my_vertices, my_bordervertices, iomodel, iomodel_handle); @@ -73,5 +77,5 @@ /*First: add all the nodes of all the elements belonging to this cpu*/ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ for (i=0;inumberofelements;i++){ if (my_elements[i]){ Index: sm/trunk/src/c/modules/ModelProcessorx/Partitioning.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Partitioning.cpp (revision 3983) +++ (revision ) @@ -1,274 +1,0 @@ -/*!\file: Partitioning.cpp - * \brief: partition elements and nodes and vertices - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include -#include "../../objects/objects.h" -#include "../../shared/shared.h" -#include "../../EnumDefinitions/EnumDefinitions.h" -#include "../../io/io.h" -#include "../../include/include.h" -#include "../MeshPartitionx/MeshPartitionx.h" -#include "../ModelProcessorx/ModelProcessorx.h" - -void DiscontinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); -void ContinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); - -void Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ - - if (iomodel->analysis_type==Prognostic2AnalysisEnum || iomodel->analysis_type==Balancedthickness2AnalysisEnum) - DiscontinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle); - else - ContinuousGalerkinPartitioning(pmy_elements, pmy_vertices, pmy_nodes, pmy_bordervertices, iomodel, iomodel_handle); -} - -void ContinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ - - /*as many nodes as there are vertices */ - - int i; - - extern int my_rank; - extern int num_procs; - - /*output: */ - bool* my_elements=NULL; - bool* my_vertices=NULL; - bool* my_nodes=NULL; - bool* my_bordervertices=NULL; - - /*intermediary: */ - int* epart=NULL; //element partitioning. - int* npart=NULL; //node partitioning. - int elements_width; //number of columns in elements (2d->3, 3d->6) - Vec bordervertices=NULL; - double* serial_bordervertices=NULL; - int el1,el2; - - /*Number of vertices per elements, needed to correctly retrieve data: */ - if(strcmp(iomodel->meshtype,"2d")==0) elements_width=3; //tria elements - else elements_width=6; //penta 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->numberofvertices,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofvertices2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); - - /*Free elements and elements2d: */ - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->elements2d); - - #else - /*In serial mode, epart is full of 0: all elements belong to cpu 0: */ - epart=(int*)xcalloc(iomodel->numberofelements,sizeof(int)); - #endif - - /*Deal with rifts, they have to be included into one partition only, not several: */ - IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo"); - - for(i=0;inumrifts;i++){ - el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing - el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - - /*Free rifts: */ - xfree((void**)&iomodel->riftinfo); - - /*Used later on: */ - my_vertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); - my_elements=(bool*)xcalloc(iomodel->numberofelements,sizeof(bool)); - - /*Start figuring out, out of the partition, which elements belong to this cpu: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - for (i=0;inumberofelements;i++){ - - /*!All elements have been partitioned above, only deal with elements for this cpu: */ - if(my_rank==epart[i]){ - - my_elements[i]=1; - - /*Now that we are here, we can also start building the list of vertices belonging to this cpu 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 vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices - will hold which vertices belong to this partition*/ - my_vertices[(int)*(iomodel->elements+elements_width*i+0)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+1)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+2)-1]=1; - - if(elements_width==6){ - my_vertices[(int)*(iomodel->elements+elements_width*i+3)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+4)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+5)-1]=1; - } - } - }//for (i=0;ielements); - - #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:*/ - bordervertices=NewVec(iomodel->numberofvertices); - - for (i=0;inumberofvertices;i++){ - if(my_vertices[i])VecSetValue(bordervertices,i,1,ADD_VALUES); - } - VecAssemblyBegin(bordervertices); - VecAssemblyEnd(bordervertices); - - VecToMPISerial(&serial_bordervertices,bordervertices); - - /*now go through serial_bordervertices, and booleanize it: */ - my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); - for(i=0;inumberofvertices;i++){ - if(serial_bordervertices[i]>1)my_bordervertices[i]=1; - } - #else - /*No border vertices: */ - my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); - #endif - - /*Deal with my_nodes: */ - my_nodes=(bool*)xmalloc(iomodel->numberofvertices*sizeof(bool)); - memcpy(my_nodes,my_vertices,iomodel->numberofvertices*sizeof(bool)); - - /*Free ressources:*/ - xfree((void**)&npart); - xfree((void**)&epart); - xfree((void**)&serial_bordervertices); - VecFree(&bordervertices); - - /*Assign output pointers:*/ - *pmy_elements=my_elements; - *pmy_vertices=my_vertices; - *pmy_nodes=my_nodes; - *pmy_bordervertices=my_bordervertices; -} - - -void DiscontinuousGalerkinPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ - - /*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. */ - - int i,j; - - /*output: */ - bool* my_elements=NULL; - bool* my_vertices=NULL; - bool* my_nodes=NULL; - bool* my_nodescontinuous=NULL; - bool* my_bordervertices=NULL; - - int i1,i2; - int cols; - double e1,e2; - int pos; - - /*First: get element and vertices partitioning from Continuous Galerkin: only the nodes are partitioned differently*/ - ContinuousGalerkinPartitioning(&my_elements,&my_vertices,&my_nodescontinuous,&my_bordervertices,iomodel,iomodel_handle); - xfree((void**)&my_nodescontinuous); - - /*Now we must build node partitioning - * - there are three nodes per element (discontinous) - * - for each element present of each partition, its three nodes will be in this partition - * - the edges require the dofs of the 2 nodes of each elements sharing the edge. - * if the 2 elements sharing the edge are on 2 different cpus, we must duplicate - * the two nodes that are not on the cpus so that the edge can access the dofs of - * all its 4 nodes - */ - - /*Allocate*/ - my_nodes=(bool*)xcalloc(3*iomodel->numberofelements,sizeof(int)); - - /*First: add all the nodes of all the elements belonging to this cpu*/ - if (strcmp(iomodel->meshtype,"2d")==0){ - for (i=0;inumberofelements;i++){ - if (my_elements[i]){ - my_nodes[3*i+0]=1; - my_nodes[3*i+1]=1; - my_nodes[3*i+2]=1; - } - } - } - else{ - ISSMERROR("not implemented yet"); - } - - /*Second: add all missing nodes*/ - - /*Get edges and elements*/ - IoModelFetchData(&iomodel->edges,&iomodel->numberofedges,&cols,iomodel_handle,"edges"); - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - if (cols!=4) ISSMERROR("field edges should have 4 columns"); - - /*!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 (internal edge) - * 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_elements[(int)e1] && !isnan(e2) && !my_elements[(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 ((int)iomodel->elements[3*(int)e2+j]==i1) pos=j; - } - - /*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){ - my_nodes[(int)e2*3+0]=1; - my_nodes[(int)e2*3+2]=1; - } - else if(pos==1){ - my_nodes[(int)e2*3+1]=1; - my_nodes[(int)e2*3+0]=1; - } - else if (pos==2){ - my_nodes[(int)e2*3+2]=1; - my_nodes[(int)e2*3+1]=1; - } - else{ - ISSMERROR("Problem in edges creation"); - } - } - } - - /*Free data: */ - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->edges); - - /*Assign output pointers:*/ - *pmy_elements=my_elements; - *pmy_vertices=my_vertices; - *pmy_nodes=my_nodes; - *pmy_bordervertices=my_bordervertices; -} Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic/CreateElementsNodesAndMaterialsPrognostic.cpp (revision 3984) @@ -34,5 +34,5 @@ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ /*Fetch data needed: */ @@ -74,5 +74,5 @@ } - else{ // if (strcmp(meshtype,"2d")==0) + else{ // if (dim==2) /*Fetch data needed: */ @@ -123,5 +123,5 @@ - } //if (strcmp(meshtype,"2d")==0) + } //if (dim==2) /*Add new constrant material property to materials, at the end: */ @@ -129,5 +129,5 @@ /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3984) @@ -37,5 +37,5 @@ /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ /*Fetch data needed: */ @@ -75,7 +75,7 @@ } - else{ // if (strcmp(meshtype,"2d")==0) + else{ // if (dim==2) ISSMERROR("not implemented yet"); - } //if (strcmp(meshtype,"2d")==0) + } //if (dim==2) /*Add new constrant material property tgo materials, at the end: */ @@ -93,5 +93,5 @@ IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateElementsNodesAndMaterialsSlopeCompute.cpp (revision 3984) @@ -34,5 +34,5 @@ /*2d mesh: */ - if (strcmp(iomodel->meshtype,"2d")==0){ + if (iomodel->dim==2){ /*Fetch data needed: */ @@ -63,5 +63,5 @@ } - else{ // if (strcmp(meshtype,"2d")==0) + else{ // if (dim==2) /*Fetch data needed: */ @@ -94,5 +94,5 @@ xfree((void**)&iomodel->lowerelements); - } //if (strcmp(meshtype,"2d")==0) + } //if (dim==2) /*Add new constrant material property tgo materials, at the end: */ @@ -100,5 +100,5 @@ /*First fetch data: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 3984) @@ -23,5 +23,5 @@ /*return if 2d mesh*/ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; + if (iomodel->dim==2)goto cleanup_and_return; /*Fetch data: */ Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateElementsNodesAndMaterialsThermal.cpp (revision 3984) @@ -96,5 +96,5 @@ /*Create nodes and vertices: */ - if (strcmp(iomodel->meshtype,"3d")==0){ + if (iomodel->dim==3){ IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); Index: /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 3983) +++ /issm/trunk/src/c/modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp (revision 3984) @@ -19,5 +19,5 @@ /*return if 2d mesh*/ - if (strcmp(iomodel->meshtype,"2d")==0)goto cleanup_and_return; + if (iomodel->dim==2)goto cleanup_and_return; /*Create loads: */ Index: sm/trunk/src/c/modules/ModelProcessorx/VerticesPartitioning.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/VerticesPartitioning.cpp (revision 3983) +++ (revision ) @@ -1,144 +1,0 @@ -/*!\file: VerticesPartitioning.cpp - * \brief: partition elements and nodes and vertices - */ - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include -#include "../../objects/objects.h" -#include "../../shared/shared.h" -#include "../../EnumDefinitions/EnumDefinitions.h" -#include "../../io/io.h" -#include "../../include/include.h" -#include "../MeshPartitionx/MeshPartitionx.h" -#include "../ModelProcessorx/ModelProcessorx.h" - -void VerticesPartitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ - - int i; - - extern int my_rank; - extern int num_procs; - - /*output: */ - bool* my_elements=NULL; - bool* my_vertices=NULL; - bool* my_bordervertices=NULL; - - /*intermediary: */ - int* epart=NULL; //element partitioning. - int* npart=NULL; //node partitioning. - int elements_width; //number of columns in elements (2d->3, 3d->6) - Vec bordervertices=NULL; - double* serial_bordervertices=NULL; - int el1,el2; - - /*Number of vertices per elements, needed to correctly retrieve data: */ - if(strcmp(iomodel->meshtype,"2d")==0) elements_width=3; //tria elements - else elements_width=6; //penta 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->numberofvertices,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofvertices2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); - - /*Free elements and elements2d: */ - xfree((void**)&iomodel->elements); - xfree((void**)&iomodel->elements2d); - - #else - /*In serial mode, epart is full of 0: all elements belong to cpu 0: */ - epart=(int*)xcalloc(iomodel->numberofelements,sizeof(int)); - #endif - - /*Deal with rifts, they have to be included into one partition only, not several: */ - IoModelFetchData(&iomodel->riftinfo,&iomodel->numrifts,NULL,iomodel_handle,"riftinfo"); - - for(i=0;inumrifts;i++){ - el1=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+2)-1; //matlab indexing to c indexing - el2=(int)*(iomodel->riftinfo+RIFTINFOSIZE*i+3)-1; //matlab indexing to c indexing - epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding grids; - } - - /*Free rifts: */ - xfree((void**)&iomodel->riftinfo); - - /*Used later on: */ - my_vertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); - my_elements=(bool*)xcalloc(iomodel->numberofelements,sizeof(bool)); - - /*Start figuring out, out of the partition, which elements belong to this cpu: */ - IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); - for (i=0;inumberofelements;i++){ - - /*!All elements have been partitioned above, only deal with elements for this cpu: */ - if(my_rank==epart[i]){ - - my_elements[i]=1; - - /*Now that we are here, we can also start building the list of vertices belonging to this cpu 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 vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices - will hold which vertices belong to this partition*/ - my_vertices[(int)*(iomodel->elements+elements_width*i+0)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+1)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+2)-1]=1; - - if(elements_width==6){ - my_vertices[(int)*(iomodel->elements+elements_width*i+3)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+4)-1]=1; - my_vertices[(int)*(iomodel->elements+elements_width*i+5)-1]=1; - } - } - }//for (i=0;ielements); - - #ifdef _PARALLEL_ - /*From the element partitioning, we can determine which vertices are on the inside of this cpu's - *element partition, and which are on its border with other vertices:*/ - bordervertices=NewVec(iomodel->numberofvertices); - - for (i=0;inumberofvertices;i++){ - if(my_vertices[i])VecSetValue(bordervertices,i,1,ADD_VALUES); - } - VecAssemblyBegin(bordervertices); - VecAssemblyEnd(bordervertices); - - VecToMPISerial(&serial_bordervertices,bordervertices); - - /*now go through serial_bordervertices, and booleanize it: */ - my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); - for(i=0;inumberofvertices;i++){ - if(serial_bordervertices[i]>1)my_bordervertices[i]=1; - } - #else - /*No border vertices: */ - my_bordervertices=(bool*)xcalloc(iomodel->numberofvertices,sizeof(bool)); - #endif - - /*Free ressources:*/ - xfree((void**)&npart); - xfree((void**)&epart); - xfree((void**)&serial_bordervertices); - VecFree(&bordervertices); - - /*Assign output pointers:*/ - *pmy_elements=my_elements; - *pmy_vertices=my_vertices; - *pmy_bordervertices=my_bordervertices; -} - - Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 3983) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 3984) @@ -27,4 +27,5 @@ this->matpar=NULL; this->neighbors=NULL; + this->collapse=NULL; this->inputs=NULL; @@ -36,15 +37,14 @@ delete inputs; this->parameters=NULL; -} -/*}}}*/ -/*FUNCTION Penta::Penta(int id, int index, IoModel* iomodel) {{{1*/ -Penta::Penta(int penta_id, int index, IoModel* iomodel){ //i is the element index + xfree((void**)&collapse); +} +/*}}}*/ +/*FUNCTION Penta::Penta(int id, int index, IoModel* iomodel,int nummodels) {{{1*/ +Penta::Penta(int penta_id, int index, IoModel* iomodel,int nummodels){ //i is the element index IssmInt i; - int penta_node_ids[6]; int penta_matice_id; int penta_matpar_id; int penta_elements_ids[2]; - double nodeinputs[6]; bool collapse; @@ -53,12 +53,7 @@ /*hooks: */ - ISSMASSERT(iomodel->elements); - for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. - penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab - } penta_matice_id=index+1; //refers to the corresponding ice material object penta_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object - ISSMASSERT(iomodel->upperelements); if isnan(iomodel->upperelements[index]){ penta_elements_ids[0]=this->id; //upper penta is the same penta @@ -67,5 +62,5 @@ penta_elements_ids[0]=(int)(iomodel->upperelements[index]); } - ISSMASSERT(iomodel->lowerelements); + if isnan(iomodel->lowerelements[index]){ penta_elements_ids[1]=this->id; //lower penta is the same penta @@ -74,8 +69,30 @@ penta_elements_ids[1]=(int)(iomodel->lowerelements[index]); } - this->InitHookNodes(penta_node_ids); this->nodes=NULL; + this->InitHookNodes(nummodels);this->nodes=NULL; this->InitHookMatice(penta_matice_id);this->matice=NULL; this->InitHookMatpar(penta_matpar_id);this->matpar=NULL; this->InitHookNeighbors(penta_elements_ids);this->neighbors=NULL; + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + + //collpase flags + collapse=(bool*)xcalloc(nummodels*sizeof(bool)); + +} +/*}}}*/ +/*FUNCTION Penta::Update(IoModel* iomodel,int analysis_counter) {{{1*/ +Penta::Update(IoModel* iomodel,int analysis_counter){ + + IssmInt i; + int penta_node_ids[6]; + double nodeinputs[6]; + + /*hooks: */ + for(i=0;i<6;i++){ //go recover node ids, needed to initialize the node hook. + penta_node_ids[i]=(int)*(iomodel->elements+6*index+i); //ids for vertices are in the elements array from Matlab + } + + this->SettHookNodes(penta_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type //intialize inputs, and add as many inputs per element as requested: @@ -186,12 +203,10 @@ if(iomodel->elements_type){ if (*(iomodel->elements_type+2*i+0)==MacAyealFormulationEnum){ - collapse=1; + collapse[analysis_counter]=true; } else{ - collapse=0; + collapse[analysis_counter]=false; } } - this->inputs->AddInput(new BoolInput(CollapseEnum,(IssmBool)collapse)); - } /*}}}*/ @@ -215,4 +230,6 @@ /*now deal with hooks and objects: */ + penta->InitHookNodes(this->nummodels); + for(i=0;inummodels;i++)penta->hnodes[i].copy(&this->hnodes[i]); penta->hnodes.copy(&this->hnodes); penta->hmatice.copy(&this->hmatice); @@ -221,5 +238,6 @@ /*recover objects: */ - penta->nodes=(Node**)penta->hnodes.deliverp(); + penta->nodes=(Node**)xmalloc(6*sizeof(Node*)); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. + for(i=0;i<6;i++)penta->nodes[i]=this->nodes[i]; penta->matice=(Matice*)penta->hmatice.delivers(); penta->matpar=(Matpar*)penta->hmatpar.delivers(); @@ -232,9 +250,9 @@ /*Object management: */ /*FUNCTION Penta::Configure {{{1*/ -void Penta::Configure(DataSet* elementsin, DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ +void Penta::Configure(DataSet* elementsin, DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin,int analysis_counter){ /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective * datasets, using internal ids and offsets hidden in hooks: */ - this->hnodes.configure(nodesin); + this->hnodes[analysis_counter].configure(nodesin); this->hmatice.configure(materialsin); this->hmatpar.configure(materialsin); @@ -242,5 +260,5 @@ /*Now, go pick up the objects inside the hooks: */ - this->nodes=(Node**)this->hnodes.deliverp(); + this->nodes=(Node**)this->hnodes[analysis_counter].deliverp(); this->matice=(Matice*)this->hmatice.delivers(); this->matpar=(Matpar*)this->hmatpar.delivers(); @@ -264,12 +282,17 @@ memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(&nummodels,marshalled_dataset,sizeof(nummodels));marshalled_dataset+=sizeof(nummodels); + + /*allocate dynamic memory: */ + collapse=(bool*)xmalloc(nummodels*sizeof(bool)); + InitHookNodes(nummodels); /*demarshall hooks: */ - hnodes.Demarshall(&marshalled_dataset); + for(i=0;iparameters=NULL; @@ -291,4 +317,6 @@ void Penta::DeepEcho(void){ + int i; + printf("Penta:\n"); printf(" id: %i\n",id); @@ -306,5 +334,6 @@ printf(" inputs\n"); inputs->DeepEcho(); - + printf(" collapse: \n "); + for(i=0;iGetParameterValue(&onsurface,ElementOnSurfaceEnum); - return onsurface; - -} -/*}}}*/ -/*FUNCTION Penta::GetUpperElement{{{1*/ -Penta* Penta::GetUpperElement(void){ - - Penta* upper_penta=NULL; - upper_penta=(Penta*)neighbors[1]; //first one under, second one above - return upper_penta; - -} -/*}}}*/ /*FUNCTION Penta::Marshall {{{1*/ void Penta::Marshall(char** pmarshalled_dataset){ + int i; char* marshalled_dataset=NULL; int enum_type=0; @@ -383,7 +372,8 @@ /*marshall Penta data: */ memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(marshalled_dataset,&nummodels,sizeof(nummodels));marshalled_dataset+=sizeof(nummodels); /*Marshall hooks: */ - hnodes.Marshall(&marshalled_dataset); + for(i=0;iMarshallSize() + +nummodels*sizeof(bool) +sizeof(int); //sizeof(int) for enum type } @@ -423,7 +422,4 @@ int indices[3]; int zero=0; - Hook* tria_hnodes=NULL; - Hook* tria_hmatice=NULL; - Hook* tria_hmatpar=NULL; Parameters* tria_parameters=NULL; Inputs* tria_inputs=NULL; @@ -433,7 +429,4 @@ indices[2]=g2; - tria_hnodes =this->hnodes.Spawn(indices,3); - tria_hmatice=this->hmatice.Spawn(&zero,1); - tria_hmatpar=this->hmatpar.Spawn(&zero,1); tria_parameters=this->parameters; tria_inputs=(Inputs*)this->inputs->SpawnTriaInputs(indices); @@ -444,17 +437,9 @@ tria->parameters=tria_parameters; - /*now deal with hooks and objects: */ - tria->hnodes.copy(tria_hnodes); - tria->hmatice.copy(tria_hmatice); - tria->hmatpar.copy(tria_hmatpar); - - /*recover objects: */ - tria->nodes=(Node**)tria->hnodes.deliverp(); - tria->matice=(Matice*)tria->hmatice.delivers(); - tria->matpar=(Matpar*)tria->hmatpar.delivers(); - - delete tria_hnodes; - delete tria_hmatice; - delete tria_hmatpar; + /*now deal with nodes, matice and matpar: */ + tria->nodes=(Node**)xmalloc(3*sizeof(Node*)); + for(i=0;i<3;i++)tria->nodes[i]=this->nodes[indices[i]]; + tria->matice=this->matice; + tria->matpar=this->matpar; return tria; @@ -468,7 +453,4 @@ int indices[2]; int zero=0; - Hook *beam_hnodes = NULL; - Hook *beam_hmatice = NULL; - Hook *beam_hmatpar = NULL; Parameters *beam_parameters = NULL; Inputs *beam_inputs = NULL; @@ -477,7 +459,4 @@ indices[1]=g1; - beam_hnodes =this->hnodes.Spawn(indices,2); - beam_hmatice=this->hmatice.Spawn(&zero,1); - beam_hmatpar=this->hmatpar.Spawn(&zero,1); beam_parameters=this->parameters; beam_inputs=(Inputs*)this->inputs->SpawnBeamInputs(indices); @@ -488,17 +467,9 @@ beam->parameters=beam_parameters; - /*now deal with hooks and objects: */ - beam->hnodes.copy(beam_hnodes); - beam->hmatice.copy(beam_hmatice); - beam->hmatpar.copy(beam_hmatpar); - - /*recover objects: */ - //beam->nodes=(Node**)beam->hnodes.deliverp(); - //beam->matice=(Matice*)beam->hmatice.delivers(); - //beam->matpar=(Matpar*)beam->hmatpar.delivers(); - - delete beam_hnodes; - delete beam_hmatice; - delete beam_hmatpar; + /*now deal with ndoes,matice and matpar: */ + beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); + for(i=0;i<2;i++)beam->nodes[i]=this->nodes[indices[i]]; + beam->matice=this->matice; + beam->matpar=this->matpar; return beam; @@ -510,13 +481,7 @@ Sing* sing=NULL; int zero=0; - Hook *sing_hnodes = NULL; - Hook *sing_hmatice = NULL; - Hook *sing_hmatpar = NULL; Parameters *sing_parameters = NULL; Inputs *sing_inputs = NULL; - sing_hnodes =this->hnodes.Spawn(&index,1); - sing_hmatice=this->hmatice.Spawn(&zero,1); - sing_hmatpar=this->hmatpar.Spawn(&zero,1); sing_parameters=this->parameters; sing_inputs=(Inputs*)this->inputs->SpawnSingInputs(index); @@ -527,18 +492,9 @@ sing->parameters=sing_parameters; - /*now deal with hooks and objects: */ - sing->hnodes.copy(sing_hnodes); - sing->hmatice.copy(sing_hmatice); - sing->hmatpar.copy(sing_hmatpar); - - /*recover objects: */ - //sing->nodes=(Node**)sing->hnodes.deliverp(); - //sing->matice=(Matice*)sing->hmatice.delivers(); - //sing->matpar=(Matpar*)sing->hmatpar.delivers(); - - delete sing_hnodes; - delete sing_hmatice; - delete sing_hmatpar; - + /*now deal with nodes,matice and matpar: */ + sing->node=this->nodes[index]; + sing->matice=this->matice; + sing->matpar=this->matpar; + return sing; } @@ -840,4 +796,45 @@ /*Object functions*/ +/*FUNCTION Penta::IsInput{{{1*/ +bool Penta::IsInput(int name){ + if (name==ThicknessEnum || + name==SurfaceEnum || + name==BedEnum || + name==SurfaceSlopexEnum || + name==SurfaceSlopeyEnum || + name==MeltingRateEnum || + name==AccumulationRateEnum || + name==GeothermalFluxEnum || + name==PressureEnum || + name==VxEnum || + name==VyEnum || + name==VzEnum || + name==TemperatureEnum || + name==TemperatureAverageEnum || + name==RheologyBEnum || + name==RheologyNEnum) { + return true; + } + else return false; +} +/*}}}*/ +/*FUNCTION Penta::IsOnSurface{{{1*/ +bool Penta::IsOnSurface(void){ + + bool onsurface; + inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum); + return onsurface; + +} +/*}}}*/ +/*FUNCTION Penta::GetUpperElement{{{1*/ +Penta* Penta::GetUpperElement(void){ + + Penta* upper_penta=NULL; + upper_penta=(Penta*)neighbors[1]; //first one under, second one above + return upper_penta; + +} +/*}}}*/ /*FUNCTION Penta::ComputeBasalStress {{{1*/ void Penta::ComputeBasalStress(Vec sigma_b,int analysis_type,int sub_analysis_type){ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 3983) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 3984) @@ -37,11 +37,15 @@ Inputs *inputs; + /*internal parameters: */ + bool *collapse; //collapse elements, depending on analysis_type + /*FUNCTION constructors, destructors {{{1*/ Penta(); - Penta(int penta_id,int i, IoModel* iomodel); + Penta(int penta_id,int i, IoModel* iomodel,int nummodels); + Update(IoModel* iomodel,int analysis_counter,int analysis_type); ~Penta(); /*}}}*/ /*FUNCTION object management {{{1*/ - void Configure(DataSet* elements,DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters); + void Configure(DataSet* elements,DataSet* loads,DataSet* nodes,DataSet* materials,Parameters* parameters,int analysis_counter); Object* copy(); void DeepEcho(); Index: /issm/trunk/src/c/objects/Elements/PentaHook.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaHook.cpp (revision 3983) +++ /issm/trunk/src/c/objects/Elements/PentaHook.cpp (revision 3984) @@ -23,13 +23,21 @@ /*FUNCTION PentaHook::PentaHook(){{{1*/ PentaHook::PentaHook(){ + numanalyses=UNDEF; + this->hooks=NULL; } /*}}}*/ /*FUNCTION PentaHook::~PentaHook(){{{1*/ PentaHook::~PentaHook(){ + xfree((void**)&hnodes); } /*}}}*/ -/*FUNCTION PentaHook::InitHookNodes(int* node_ids){{{1*/ -void PentaHook::InitHookNodes(int* node_ids){ - this->hnodes.Init(node_ids,6); +/*FUNCTION PentaHook::InitHookNodes(int numanalyses){{{1*/ +void PentaHook::InitHookNodes(int numanalyses){ + this->hnodes=(Hook*)xmalloc(numanalyses*sizeof(Hook)); +} +/*}}}*/ +/*FUNCTION PentaHook::SetHookNodes(int* node_ids,int analysis_counter){{{1*/ +void PentaHook::SetHookNodes(int* node_ids,int analysis_counter){ + this->hnodes[analysis_counter].Init(node_ids,6); } Index: /issm/trunk/src/c/objects/Elements/PentaHook.h =================================================================== --- /issm/trunk/src/c/objects/Elements/PentaHook.h (revision 3983) +++ /issm/trunk/src/c/objects/Elements/PentaHook.h (revision 3984) @@ -11,5 +11,6 @@ public: - Hook hnodes; // 6 nodes + int numanalyses; //number of analysis types + Hook* hnodes; // 6 nodes for each analysis type Hook hmatice; // 1 ice material Hook hmatpar; // 1 material parameter @@ -20,5 +21,6 @@ PentaHook(); ~PentaHook(); - void InitHookNodes(int* node_ids); + void InitHookNodes(int numanalyses); + void SetHookNodes(int* node_ids,int analysis_counter); void InitHookMatice(int matice_id); void InitHookMatpar(int matpar_id); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3983) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 3984) @@ -28,4 +28,5 @@ this->inputs=NULL; this->parameters=NULL; + this->collapse=NULL; } /*}}}*/ @@ -34,15 +35,13 @@ delete inputs; this->parameters=NULL; -} -/*}}}*/ -/*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel){{{1*/ -Tria::Tria(int tria_id, int index, IoModel* iomodel){ //i is the element index + xfree((void**)&collapse); +} +/*}}}*/ +/*FUNCTION Tria::Tria(int id, int index, IoModel* iomodel,int nummodels){{{1*/ +Tria::Tria(int tria_id, int index, IoModel* iomodel,int nummodels){ //i is the element index int i; - int j; - int tria_node_ids[3]; int tria_matice_id; int tria_matpar_id; - double nodeinputs[3]; /*id: */ @@ -52,4 +51,25 @@ this->InitInterpolationType(P1Enum); + /*hooks: */ + tria_matice_id=index+1; //refers to the corresponding ice material object + tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object + + this->InitHookNodes(nummodels);this->nodes=NULL; + this->InitHookNodes(tria_node_ids); this->nodes=NULL; + this->InitHookMatice(tria_matice_id);this->matice=NULL; + this->InitHookMatpar(tria_matpar_id);this->matpar=NULL; + + //this->parameters: we still can't point to it, it may not even exist. Configure will handle this. + this->parameters=NULL; + +} +/*}}}*/ +/*FUNCTION Tria::Update(IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/ +Tria::Update(IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index + + int i; + int tria_node_ids[3]; + double nodeinputs[3]; + /*hooks: */ //go recover node ids, needed to initialize the node hook. @@ -67,10 +87,6 @@ } } - tria_matice_id=index+1; //refers to the corresponding ice material object - tria_matpar_id=iomodel->numberofelements+1; //refers to the constant material parameters object - - this->InitHookNodes(tria_node_ids); this->nodes=NULL; - this->InitHookMatice(tria_matice_id);this->matice=NULL; - this->InitHookMatpar(tria_matpar_id);this->matpar=NULL; + + this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type //intialize inputs, and add as many inputs per element as requested: @@ -204,4 +220,6 @@ /*now deal with hooks and objects: */ + tria->InitHookNodes(this->nummodels); + for(i=0;inummodels;i++)tria->hnodes[i].copy(&this->hnodes[i]); tria->hnodes.copy(&this->hnodes); tria->hmatice.copy(&this->hmatice); @@ -209,4 +227,6 @@ /*recover objects: */ + tria->nodes=(Node**)xmalloc(3*sizeof(Node*)); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes. + for(i=0;i<3;i++)tria->nodes[i]=this->nodes[i]; tria->nodes=(Node**)tria->hnodes.deliverp(); tria->matice=(Matice*)tria->hmatice.delivers(); @@ -219,14 +239,14 @@ /*Object management: */ /*FUNCTION Tria::Configure {{{1*/ -void Tria::Configure(DataSet* elementsin, DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin){ +void Tria::Configure(DataSet* elementsin, DataSet* loadsin, DataSet* nodesin, DataSet* materialsin, Parameters* parametersin,int analysis_counter){ /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective * datasets, using internal ids and offsets hidden in hooks: */ - this->hnodes.configure(nodesin); + this->hnodes[analysis_counter].configure(nodesin); this->hmatice.configure(materialsin); this->hmatpar.configure(materialsin); /*Now, go pick up the objects inside the hooks: */ - this->nodes=(Node**)this->hnodes.deliverp(); + this->nodes=(Node**)this->hnodes[analysis_counter].deliverp(); this->matice=(Matice*)this->hmatice.delivers(); this->matpar=(Matpar*)this->hmatpar.delivers(); @@ -241,4 +261,5 @@ char* marshalled_dataset=NULL; + int i; /*recover marshalled_dataset: */ @@ -249,7 +270,11 @@ memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); memcpy(&interpolation_type,marshalled_dataset,sizeof(interpolation_type));marshalled_dataset+=sizeof(interpolation_type); + memcpy(&nummodels,marshalled_dataset,sizeof(nummodels));marshalled_dataset+=sizeof(nummodels); + + /*allocate dynamic memory: */ + InitHookNodes(nummodels); /*demarshall hooks: */ - hnodes.Demarshall(&marshalled_dataset); + for(i=0;ihnodes.Spawn(indices,2); - beam_hmatice=this->hmatice.Spawn(&zero,1); - beam_hmatpar=this->hmatpar.Spawn(&zero,1); beam_parameters=this->parameters; beam_inputs=(Inputs*)this->inputs->SpawnBeamInputs(indices); @@ -395,17 +413,10 @@ beam->parameters=beam_parameters; - /*now deal with hooks and objects: */ - beam->hnodes.copy(beam_hnodes); - beam->hmatice.copy(beam_hmatice); - beam->hmatpar.copy(beam_hmatpar); - - /*recover objects: */ - //beam->nodes=(Node**)beam->hnodes.deliverp(); - //beam->matice=(Matice*)beam->hmatice.delivers(); - //beam->matpar=(Matpar*)beam->hmatpar.delivers(); - - delete beam_hnodes; - delete beam_hmatice; - delete beam_hmatpar; + /*now deal with nodes, matice and matpar: */ + beam->nodes=(Node**)xmalloc(2*sizeof(Node*)); + for(i=0;i<2;i++)beam->nodes[i]=this->nodes[indices[i]]; + beam->matice=this->matice; + beam->matpar=this->matpar; + return beam; @@ -417,13 +428,7 @@ Sing* sing=NULL; int zero=0; - Hook *sing_hnodes = NULL; - Hook *sing_hmatice = NULL; - Hook *sing_hmatpar = NULL; Parameters *sing_parameters = NULL; Inputs *sing_inputs = NULL; - sing_hnodes =this->hnodes.Spawn(&index,1); - sing_hmatice=this->hmatice.Spawn(&zero,1); - sing_hmatpar=this->hmatpar.Spawn(&zero,1); sing_parameters=this->parameters; sing_inputs=(Inputs*)this->inputs->SpawnSingInputs(index); @@ -434,18 +439,9 @@ sing->parameters=sing_parameters; - /*now deal with hooks and objects: */ - sing->hnodes.copy(sing_hnodes); - sing->hmatice.copy(sing_hmatice); - sing->hmatpar.copy(sing_hmatpar); - - /*recover objects: */ - //sing->nodes=(Node**)sing->hnodes.deliverp(); - //sing->matice=(Matice*)sing->hmatice.delivers(); - //sing->matpar=(Matpar*)sing->hmatpar.delivers(); - - delete sing_hnodes; - delete sing_hmatice; - delete sing_hmatpar; - + /*now deal with node,matice and matpar: */ + sing->node=this->nodes[index]; + sing->matice=this->matice; + sing->matpar=this->matpar; + return sing; } @@ -686,4 +682,13 @@ /*Object functions*/ +/*FUNCTION Tria::IsInput{{{1*/ +bool Tria::IsInput(int name){ + if (name==SurfaceSlopexEnum || + name==SurfaceSlopeyEnum){ + return true; + } + else return false; +} +/*}}}*/ /*FUNCTION Tria::ComputeBasalStress {{{1*/ void Tria::ComputeBasalStress(Vec eps,int analysis_type,int sub_analysis_type){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 3983) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 3984) @@ -38,4 +38,5 @@ Tria(); Tria(int tria_id,int i, IoModel* iomodel); + Update(IoModel* iomodel,int analysis_counter,int analysis_type); ~Tria(); /*}}}*/ Index: /issm/trunk/src/c/objects/Elements/TriaHook.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaHook.cpp (revision 3983) +++ /issm/trunk/src/c/objects/Elements/TriaHook.cpp (revision 3984) @@ -23,13 +23,21 @@ /*FUNCTION TriaHook::TriaHook(){{{1*/ TriaHook::TriaHook(){ + numanalyses=UNDEF; + this->hooks=NULL; } /*}}}*/ /*FUNCTION TriaHook::~TriaHook(){{{1*/ TriaHook::~TriaHook(){ + xfree((void**)&hnodes); } /*}}}*/ -/*FUNCTION TriaHook::InitHookNodes(int* node_ids){{{1*/ -void TriaHook::InitHookNodes(int* node_ids){ - this->hnodes.Init(node_ids,3); +/*FU/*FUNCTION TriaHook::InitHookNodes(int numanalyses){{{1*/ +void TriaHook::InitHookNodes(int numanalyses){ + this->hnodes=(Hook*)xmalloc(numanalyses*sizeof(Hook)); +} +/*}}}*/ +/*FUNCTION TriaHook::SetHookNodes(int* node_ids,int analysis_counter){{{1*/ +void TriaHook::SetHookNodes(int* node_ids,int analysis_counter){ + this->hnodes[analysis_counter].Init(node_ids,3); } Index: /issm/trunk/src/c/objects/Elements/TriaHook.h =================================================================== --- /issm/trunk/src/c/objects/Elements/TriaHook.h (revision 3983) +++ /issm/trunk/src/c/objects/Elements/TriaHook.h (revision 3984) @@ -11,5 +11,6 @@ public: - Hook hnodes; // 3 nodes + int numanalyses; //number of analysis types + Hook* hnodes; // 3 nodes for each analysis type Hook hmatice; // 1 ice material Hook hmatpar; // 1 material parameter @@ -19,7 +20,10 @@ TriaHook(); ~TriaHook(); - void InitHookNodes(int* node_ids); + void InitHookNodes(int numanalyses); + void SetHookNodes(int* node_ids,int analysis_counter); void InitHookMatice(int matice_id); void InitHookMatpar(int matpar_id); + void InitHookNeighbors(int* element_ids); + /*}}}*/ Index: /issm/trunk/src/c/objects/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 3983) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 3984) @@ -24,5 +24,5 @@ nummodels=in_nummodels; - current_analysis_type=-1; + analysis_counter=-1; /*Dynamically allocate whatever is a list of length nummodels: */ @@ -97,5 +97,5 @@ for(i=0;ivector,nodesets[i]); + Reducevectorgtosx(&ys[analysis_counter], yg[analysis_counter]->vector,nodesets[analysis_counter]); _printf_(" normalizing rigid body constraints matrix:\n"); - NormalizeConstraintsx(&Gmn[i], Rmg[i],nodesets[i]); + NormalizeConstraintsx(&Gmn[analysis_counter], Rmg[analysis_counter],nodesets[analysis_counter]); _printf_(" configuring element and loads:\n"); @@ -155,5 +150,5 @@ /*FUNCTION FemModel::GetCurrentAnalysis {{{1*/ FemFemModel* FemModel::GetCurrentAnalysis(){ - return analysis_type_list[current_analysis_type]; + return analysis_type_list[analysis_counter]; } /*}}}1*/ @@ -167,5 +162,5 @@ } } - if(found)current_analysis_type=i; + if(found)analysis_counter=i; else ISSMERRR("%s%s%s"," could not find analysis_type ",EnumAsString(analysis_type) " in list of FemModel analyses"); } Index: /issm/trunk/src/c/objects/FemModel.h =================================================================== --- /issm/trunk/src/c/objects/FemModel.h (revision 3983) +++ /issm/trunk/src/c/objects/FemModel.h (revision 3984) @@ -24,6 +24,6 @@ int nummodels; - int* analysis_type_list; //list of analyses this femmodel is capable of. - int current_analysis_type; //counter to the current analysis being carried out + int* analysis_type_list; //list of analyses this femmodel is going to carry out + int analysis_counter; //counter into analysis_type_list DataSet* elements; //elements (one set for all analyses) @@ -38,5 +38,5 @@ DofVec* tpartition; - //multiple sets of vectors for each analysis_type + //multiple sets of matrices/vectors for each analysis_type Mat* Rmg; //rigid body motions matrices Mat* Gmn; Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 3983) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 3984) @@ -41,5 +41,5 @@ xfree((void**)&this->gridonhutter); xfree((void**)&this->gridonmacayeal); - if (strcmp(this->meshtype,"3d")==0){ + if (this->dim==3){ xfree((void**)&this->elements2d); xfree((void**)&this->deadgrids); @@ -96,5 +96,4 @@ xfree((void**)&this->outputfilename); xfree((void**)&this->repository); - xfree((void**)&this->meshtype); xfree((void**)&this->name); @@ -132,10 +131,10 @@ IoModelFetchData(&this->qmu_analysis,iomodel_handle,"qmu_analysis"); IoModelFetchData(&this->control_analysis,iomodel_handle,"control_analysis"); - IoModelFetchData(&this->meshtype,iomodel_handle,"type"); + IoModelFetchData(&this->dim,iomodel_handle,"dim"); /*!Get numberofelements and numberofvertices: */ IoModelFetchData(&this->numberofvertices,iomodel_handle,"numberofgrids"); IoModelFetchData(&this->numberofelements,iomodel_handle,"numberofelements"); /*!In case we are running 3d, we are going to need the collapsed and non-collapsed 2d meshes, from which the 3d mesh was extruded: */ - if (strcmp(this->meshtype,"3d")==0){ + if (this->dim==3){ /*!Deal with 2d mesh: */ @@ -229,5 +228,4 @@ this->outputfilename=NULL; this->repository=NULL; - this->meshtype=NULL; this->qmu_analysis=0; this->control_analysis=0; @@ -392,5 +390,5 @@ int i,j; - if(which_part==1 && my_rank==rank && (strcmp(this->meshtype,"3d")==0)){ + if(which_part==1 && my_rank==rank && this->dim==3){ printf("IoModel penalties: \n"); printf(" number of penalties: %i\n",this->numpenalties); Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 3983) +++ /issm/trunk/src/c/objects/IoModel.h (revision 3984) @@ -19,5 +19,5 @@ char* outputfilename; char* repository; - char* meshtype; + int dim; int qmu_analysis; int control_analysis; Index: /issm/trunk/src/c/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 3983) +++ /issm/trunk/src/c/objects/Materials/Matice.cpp (revision 3984) @@ -35,5 +35,5 @@ /*}}}*/ /*FUNCTION Matice::Matice(int id, int i, IoModel* iomodel, int num_vertices){{{1*/ -Matice::Matice(int matice_mid,int i, IoModel* iomodel, int num_vertices){ +Matice::Matice(int matice_mid,int i, IoModel* iomodel){ int j; @@ -45,35 +45,27 @@ /*intermediary: */ double B_avg; - + + /*2d or 3d? */ + if(strcmp(iomodel->meshtype,"2d")==0){ + num_vertices=3; //tria elements + } + else if(strcmp(iomodel->meshtype,"3d")==0){ + num_vertices=6; //penta elements + } + else ISSMERROR(" Mesh type not supported yet!"); + /*Compute B on the element if provided*/ - if (num_vertices==3 || num_vertices==6){ - if (iomodel->rheology_B){ - B_avg=0; - for(j=0;jrheology_B+((int)*(iomodel->elements+num_vertices*i+j)-1)); - } - B_avg=B_avg/num_vertices; - } - } - else if (num_vertices==1 || num_vertices==2){ - if (iomodel->rheology_B){ - B_avg=*(iomodel->rheology_B+i); - } - } - else ISSMERROR("num_vertices = %i not supported yet",num_vertices); - - if (iomodel->rheology_B) matice_B=B_avg; - else matice_B=UNDEF; - if (iomodel->rheology_n){ - if (num_vertices==3 || num_vertices==6){ - matice_n=(double)*(iomodel->rheology_n+i); - } - else if (num_vertices==1 || num_vertices==2){ - /*n is on the elements for now, so just take the first one*/ - matice_n=(double)*(iomodel->rheology_n); - } - else ISSMERROR("num_vertices = %i not supported yet",num_vertices); - } - else matice_n=UNDEF; + if (iomodel->rheology_B){ + B_avg=0; + for(j=0;jrheology_B+((int)*(iomodel->elements+num_vertices*i+j)-1)); + } + B_avg=B_avg/num_vertices; + matice_B=B_avg; + } + else matice_B=UNDEF; + + if (iomodel->rheology_n) matice_n=(double)*(iomodel->rheology_n+i); + else matice_n=UNDEF; this->Init(matice_mid,matice_B,matice_n); Index: /issm/trunk/src/c/objects/Materials/Matice.h =================================================================== --- /issm/trunk/src/c/objects/Materials/Matice.h (revision 3983) +++ /issm/trunk/src/c/objects/Materials/Matice.h (revision 3984) @@ -23,5 +23,5 @@ Matice(); Matice(int mid,double B,double n); - Matice(int mid,int i, IoModel* iomodel, int num_vertices); + Matice(int mid,int i, IoModel* iomodel); void Init(int mid,double B,double n); ~Matice(); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 3983) +++ /issm/trunk/src/c/objects/Node.cpp (revision 3984) @@ -36,6 +36,6 @@ } /*}}}*/ -/*FUNCTION Node::Node(int id, DofIndexing* indexing, Hook* vertex, Hook* uppernode,Inputs* inputs){{{2*/ -Node::Node(int node_id,DofIndexing* node_indexing, Hook* node_vertex, Hook* node_upper_node,Inputs* node_inputs): +/*FUNCTION Node::Node(int id, DofIndexing* indexing, Hook* vertex, Hook* uppernode,Inputs* inputs,int analysis_type){{{2*/ +Node::Node(int node_id,DofIndexing* node_indexing, Hook* node_vertex, Hook* node_upper_node,Inputs* node_inputs,int analysis_type): indexing(node_indexing), hvertex(node_vertex), @@ -44,4 +44,5 @@ /*all the initialization has been done by the initializer, just fill in the id: */ this->id=node_id; + this->analysis_type=analysis_type; if(node_inputs){ @@ -53,6 +54,6 @@ } /*}}}*/ -/*FUNCTION Node::Node(int id, int i, IoModel* iomodel) -> Continuous Galerkin{{{2*/ -Node::Node(int node_id, int i, IoModel* iomodel){ //i is the node index +/*FUNCTION Node::Node(int id, int i, IoModel* iomodel,int analysis_type) -> Continuous Galerkin{{{2*/ +Node::Node(int node_id, int i, IoModel* iomodel,int analysis_type){ //i is the node index int k; @@ -63,4 +64,5 @@ /*id: */ this->id=node_id; //matlab indexing + this->analysis_type=analysis_type; /*indexing:*/ @@ -167,6 +169,6 @@ } /*}}}*/ -/*FUNCTION Node::Node(int id, int i, int j, IoModel* iomodel) -> Discontinuous Galerkin{{{2*/ -Node::Node(int node_id,int i,int j,IoModel* iomodel){ +/*FUNCTION Node::Node(int id, int i, int j, IoModel* iomodel,int analysis_type) -> Discontinuous Galerkin{{{2*/ +Node::Node(int node_id,int i,int j,IoModel* iomodel,int analysis_type){ /* i -> index of the vertex in C indexing * j -> index of the node in C indexing*/ @@ -178,4 +180,5 @@ /*id: */ this->id=node_id; //matlab indexing + this->analysis_type=analysis_type; /*indexing:*/ @@ -236,5 +239,5 @@ Object* Node::copy() { - return new Node(this->id,&this->indexing, &this->hvertex,&this->hupper_node,this->inputs); + return new Node(this->id,&this->indexing, &this->hvertex,&this->hupper_node,this->inputs,this->analysis_type); } @@ -246,4 +249,5 @@ printf("Node:\n"); printf(" id: %i\n",id); + printf(" analysis_type: %s\n",EnumAsString(analysis_type)); indexing.DeepEcho(); printf("Vertex:\n"); @@ -269,4 +273,5 @@ memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(&analysis_type,marshalled_dataset,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type); /*demarshall objects: */ @@ -288,4 +293,5 @@ printf("Node:\n"); printf(" id: %i\n",id); + printf(" analysis_type: %s\n",EnumAsString(analysis_type)); indexing.Echo(); hvertex.Echo(); @@ -372,4 +378,5 @@ /*marshall Node data: */ memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id); + memcpy(marshalled_dataset,&analysis_type,sizeof(analysis_type));marshalled_dataset+=sizeof(analysis_type); /*marshall objects: */ @@ -399,4 +406,5 @@ hupper_node.MarshallSize()+ inputs->MarshallSize()+ + sizeof(analysis_type)+ sizeof(int); //sizeof(int) for enum type } Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 3983) +++ /issm/trunk/src/c/objects/Node.h (revision 3984) @@ -27,5 +27,6 @@ Hook hvertex; Hook hupper_node; - Inputs* inputs; //properties of this node + Inputs* inputs; //properties of this node + int analysis_type; public: @@ -34,7 +35,7 @@ Node(); Node(int id,int vertex_id, int upper_node_id, int numberofdofs); - Node(int id,DofIndexing* indexing, Hook* vertex, Hook* upper_node, Inputs* inputs); - Node(int id, int i, IoModel* iomodel); - Node(int id, int i,int j,IoModel* iomodel); + Node(int id,DofIndexing* indexing, Hook* vertex, Hook* upper_node, Inputs* inputs,int analysis_type); + Node(int id, int i, IoModel* iomodel,int analysis_type); + Node(int id, int i,int j,IoModel* iomodel,int analysis_type); ~Node(); /*}}}*/