Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3981) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3982) @@ -26,4 +26,8 @@ //diagnostic DiagnosticAnalysisEnum, + DiagnosticHorizAnalysisEnum, + DiagnosticVertAnalysisEnum, + DiagnosticHutterAnalysisEnum, + DiagnosticStokesAnalysisEnum, HorizAnalysisEnum, HutterAnalysisEnum, Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3981) +++ /issm/trunk/src/c/Makefile.am (revision 3982) @@ -78,6 +78,4 @@ ./objects/Bamg/VertexOnVertex.cpp\ ./objects/Element.h\ - ./objects/Model.h\ - ./objects/Model.cpp\ ./objects/FemModel.h\ ./objects/FemModel.cpp\ @@ -533,6 +531,4 @@ ./objects/Bamg/VertexOnVertex.h\ ./objects/Element.h\ - ./objects/Model.h\ - ./objects/Model.cpp\ ./objects/FemModel.h\ ./objects/FemModel.cpp\ Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 3981) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 3982) @@ -16,39 +16,39 @@ -void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ +void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pvertices, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,Parameters** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_type){ + + /*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); /*create parameters common to all solutions: */ CreateParameters(pparameters,iomodel,iomodel_handle); + CreateParametersControl(pparameters,iomodel,iomodel_handle); /*This is just a high level driver: */ - if (iomodel->analysis_type==DiagnosticAnalysisEnum){ + 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){ - if (iomodel->sub_analysis_type==HorizAnalysisEnum){ + CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); + CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle); + } + else if (analysis_type==DiagnosticStokesAnalysisEnum){ - CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticHoriz(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticHoriz(ploads,iomodel,iomodel_handle); - } - else if (iomodel->sub_analysis_type==VertAnalysisEnum){ - - CreateElementsNodesAndMaterialsDiagnosticVert(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticVert(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticVert(ploads,iomodel,iomodel_handle); - } - else if (iomodel->sub_analysis_type==StokesAnalysisEnum){ + CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); + CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle); + } + else if (analysis_type==DiagnosticHutterAnalysisEnum){ - CreateElementsNodesAndMaterialsDiagnosticStokes(pelements,pnodes, pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticStokes(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticStokes(ploads,iomodel,iomodel_handle); - } - else if (iomodel->sub_analysis_type==HutterAnalysisEnum){ - - CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); - CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle); - CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle); - } + CreateElementsNodesAndMaterialsDiagnosticHutter(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); + CreateConstraintsDiagnosticHutter(pconstraints,iomodel,iomodel_handle); + CreateLoadsDiagnosticHutter(ploads,iomodel,iomodel_handle); } - else if (iomodel->analysis_type==SlopecomputeAnalysisEnum){ + else if (analysis_type==SlopecomputeAnalysisEnum){ CreateElementsNodesAndMaterialsSlopeCompute(pelements,pnodes, pvertices,pmaterials, iomodel,iomodel_handle); @@ -56,18 +56,18 @@ CreateLoadsSlopeCompute(ploads,iomodel,iomodel_handle); } - else if (iomodel->analysis_type==ThermalAnalysisEnum){ + 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 (iomodel->analysis_type==MeltingAnalysisEnum){ - + 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 (iomodel->analysis_type==PrognosticAnalysisEnum){ + else if (analysis_type==PrognosticAnalysisEnum){ CreateElementsNodesAndMaterialsPrognostic(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); @@ -75,5 +75,5 @@ CreateLoadsPrognostic(ploads,iomodel,iomodel_handle); } - else if (iomodel->analysis_type==Prognostic2AnalysisEnum){ + else if (analysis_type==Prognostic2AnalysisEnum){ CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); @@ -81,5 +81,5 @@ CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle); } - else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum){ + else if (analysis_type==BalancedthicknessAnalysisEnum){ CreateElementsNodesAndMaterialsBalancedthickness(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); @@ -87,5 +87,5 @@ CreateLoadsBalancedthickness(ploads,iomodel,iomodel_handle); } - else if (iomodel->analysis_type==Balancedthickness2AnalysisEnum){ + else if (analysis_type==Balancedthickness2AnalysisEnum){ CreateElementsNodesAndMaterialsBalancedthickness2(pelements,pnodes,pvertices, pmaterials, iomodel,iomodel_handle); @@ -93,5 +93,5 @@ CreateLoadsBalancedthickness2(ploads,iomodel,iomodel_handle); } - else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum){ + else if (analysis_type==BalancedvelocitiesAnalysisEnum){ CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pvertices,pmaterials, iomodel,iomodel_handle); @@ -100,6 +100,5 @@ } else{ - - ISSMERROR("%s%i%s%i%s"," analysis_type: ",iomodel->analysis_type," sub_analysis_type: ",iomodel->sub_analysis_type," not supported yet!"); + ISSMERROR("%s%i%s%i%s"," analysis_type: ",analysis_type," analysis_type: ",analysis_type," not supported yet!"); } Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 3981) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 3982) @@ -15,5 +15,5 @@ /*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); +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); @@ -87,5 +87,6 @@ /*partitioning: */ -void Partitioning(bool** pmy_elements, bool** pmy_vertices, bool** pmy_nodes, bool** pmy_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); +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); /*Conectivity*/ Index: /issm/trunk/src/c/modules/ModelProcessorx/NodesPartitioning.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/NodesPartitioning.cpp (revision 3982) +++ /issm/trunk/src/c/modules/ModelProcessorx/NodesPartitioning.cpp (revision 3982) @@ -0,0 +1,145 @@ +/*!\file: NodesPartitioning.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 DiscontinuousGalerkinNodesPartitioning(bool* pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle); +void ContinuousGalerkinNodesPartitioning(bool** pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_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){ + + if(continuous==true){ + ContinuousGalerkinNodesPartitioning(pmy_nodes,my_elements, my_vertices, my_bordervertices, iomodel, iomodel_handle); + else + DiscontinuousGalerkinNodesPartitioning(pmy_nodes,my_elements, my_vertices, my_bordervertices, iomodel, iomodel_handle); +} + +void ContinuousGalerkinNodesPartitioning(bool** pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_bordervertices, IoModel* iomodel, ConstDataHandle iomodel_handle){ + + /*as many nodes as there are vertices */ + + /*output: */ + bool* my_nodes=NULL; + + my_nodes=(bool*)xmalloc(iomodel->numberofvertices*sizeof(bool)); + memcpy(my_nodes,my_vertices,iomodel->numberofvertices*sizeof(bool)); + + /*Assign output pointers:*/ + *pmy_nodes=my_nodes; +} + + +void DiscontinuousGalerkinNodesPartitioning(bool* pmy_nodes,bool* my_elements, bool* my_vertices, bool* my_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_nodes=NULL; + + int i1,i2; + int cols; + double e1,e2; + int pos; + + /*Build discontinuous 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_nodes=my_nodes; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/VerticesPartitioning.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/VerticesPartitioning.cpp (revision 3982) +++ /issm/trunk/src/c/modules/ModelProcessorx/VerticesPartitioning.cpp (revision 3982) @@ -0,0 +1,144 @@ +/*!\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/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 3981) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 3982) @@ -2,4 +2,5 @@ * \brief: implementation of the FemModel object */ + #ifdef HAVE_CONFIG_H @@ -10,76 +11,163 @@ #include "stdio.h" -#include "../shared/shared.h" +#include "../DataSet/DataSet.h" +#include "../modules/FemModelProcessorx/FemModelProcessorx.h" +#include "./objects.h" #include "../include/include.h" -#include "./objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" -/*constructors/destructors*/ -/*FUNCTION FemModel::FemModel {{{1*/ -FemModel::FemModel(){ - elements=NULL; - nodes=NULL; - vertices=NULL; - constraints=NULL; - loads=NULL; - materials=NULL; - parameters=NULL; +/*Object constructors and destructor*/ +/*FUNCTION FemModel::constructor {{{1*/ +FemModel::FemModel(int in_nummodels){ - partition=NULL; - tpartition=NULL; - yg=NULL; - Rmg=NULL; - nodesets=NULL; - ys=NULL; - ys0=NULL; - Gmn=NULL; + nummodels=in_nummodels; + current_analysis_type=-1; + + /*Dynamically allocate whatever is a list of length nummodels: */ + analysis_type_list=(int*)xmalloc(nummodels*sizeof(int)); + for(i=0;ivector,nodesets[i]); + + _printf_(" normalizing rigid body constraints matrix:\n"); + NormalizeConstraintsx(&Gmn[i], Rmg[i],nodesets[i]); + + _printf_(" configuring element and loads:\n"); + ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); + + _printf_(" process parameters:\n"); + ProcessParamsx( parameters, partition->vector); + + _printf_(" free ressources:\n"); + delete iomodel; +} +/*}}}1*/ +/*FUNCTION FemModel::GetCurrentAnalysis {{{1*/ +FemFemModel* FemModel::GetCurrentAnalysis(){ + return analysis_type_list[current_analysis_type]; +} +/*}}}1*/ +/*FUNCTION FemModel::SetCurrentAnalysis {{{1*/ +void FemModel::SetCurrentAnalysis(int analysis_type){ + int found=-1; + for(i=0;iinputfilename,iomodel_handle,"inputfilename"); IoModelFetchData(&this->outputfilename,iomodel_handle,"outputfilename"); - IoModelFetchData(&this->analysis_type,iomodel_handle,"analysis_type"); - IoModelFetchData(&this->sub_analysis_type,iomodel_handle,"sub_analysis_type"); IoModelFetchData(&this->qmu_analysis,iomodel_handle,"qmu_analysis"); IoModelFetchData(&this->control_analysis,iomodel_handle,"control_analysis"); @@ -232,6 +230,4 @@ this->repository=NULL; this->meshtype=NULL; - this->analysis_type=0; - this->sub_analysis_type=0; this->qmu_analysis=0; this->control_analysis=0; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 3981) +++ /issm/trunk/src/c/objects/IoModel.h (revision 3982) @@ -20,6 +20,4 @@ char* repository; char* meshtype; - int analysis_type; - int sub_analysis_type; int qmu_analysis; int control_analysis; Index: sm/trunk/src/c/objects/Model.cpp =================================================================== --- /issm/trunk/src/c/objects/Model.cpp (revision 3981) +++ (revision ) @@ -1,645 +1,0 @@ -/*!\file Model.c - * \brief: implementation of the Model object - */ - - -#ifdef HAVE_CONFIG_H - #include "config.h" -#else -#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" -#endif - -#include "stdio.h" -#include "../DataSet/DataSet.h" -#include "../modules/ModelProcessorx/ModelProcessorx.h" -#include "./objects.h" -#include "../include/include.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../modules/modules.h" - - -/*Object constructors and destructor*/ -/*FUNCTION Model::constructor {{{1*/ -Model::Model(){ - - femmodels=new DataSet(0); - active=NULL; - -} -/*}}}1*/ -/*FUNCTION Model::destructor {{{1*/ -Model::~Model(){ - - delete femmodels; - /*do not delete active, it just points to one of the femmodels!*/ -} -/*}}}1*/ - -/*Object functions*/ -/*FUNCTION Model::AddFormulation(ConstDataHandle MODEL, int analysis_type,int sub_analysis_type){{{1*/ -void Model::AddFormulation(ConstDataHandle MODEL, int analysis_type,int sub_analysis_type){ - - /*FemModel: */ - FemModel* femmodel=NULL; - DataSet* elements=NULL; - DataSet* nodes=NULL; - DataSet* vertices=NULL; - DataSet* constraints=NULL; - DataSet* loads=NULL; - DataSet* materials=NULL; - Parameters* parameters=NULL; - DofVec* partition=NULL; - DofVec* tpartition=NULL; - DofVec* yg=NULL; - Mat Rmg=NULL; - Mat Gmn=NULL; - NodeSets* nodesets=NULL; - Vec ys=NULL; - Vec ys0=NULL; - - /*intermediary: */ - IoModel* iomodel=NULL; - - _printf_(" fill model with matlab workspace data\n"); - iomodel=new IoModel(MODEL); - - _printf_(" specifying analysis\n"); - if (analysis_type!=0){ - iomodel->analysis_type=analysis_type; - } - if (sub_analysis_type!=0){ - iomodel->sub_analysis_type=sub_analysis_type; - } - - _printf_(" create datasets:\n"); - CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,MODEL); - - _printf_(" create degrees of freedom: \n"); - Dofx( &partition,&tpartition,elements,nodes, vertices,parameters); - - _printf_(" create single point constraints: \n"); - SpcNodesx( &yg, nodes,constraints); - - _printf_(" create rigid body constraints:\n"); - MpcNodesx( &Rmg, nodes,constraints); - - _printf_(" create node sets:\n"); - BuildNodeSetsx(&nodesets, nodes); - - _printf_(" reducing single point constraints vector:\n"); - Reducevectorgtosx(&ys,yg->vector,nodesets); - - _printf_(" normalizing rigid body constraints matrix:\n"); - NormalizeConstraintsx(&Gmn, Rmg,nodesets); - - _printf_(" configuring element and loads:\n"); - ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); - - _printf_(" process parameters:\n"); - ProcessParamsx( parameters, partition->vector); - - _printf_(" free ressources:\n"); - delete iomodel; - - /*Use all the data created to create an femmodel: */ - femmodel=new FemModel(elements, nodes, vertices, constraints,loads,materials,parameters, - partition,tpartition,yg,Rmg,Gmn,nodesets,ys,ys0); - - /*Add to femmodels in model: */ - femmodels->AddObject((Object*)femmodel); - -} -/*}}}1*/ -/*FUNCTION Model::AddFormulation(ConstDataHandle MODEL, int analysis_type) {{{1*/ -void Model::AddFormulation(ConstDataHandle MODEL, int analysis_type){ - - /*FemModel: */ - FemModel* femmodel=NULL; - DataSet* elements=NULL; - DataSet* nodes=NULL; - DataSet* vertices=NULL; - DataSet* constraints=NULL; - DataSet* loads=NULL; - DataSet* materials=NULL; - Parameters* parameters=NULL; - DofVec* partition=NULL; - DofVec* tpartition=NULL; - DofVec* yg=NULL; - Mat Rmg=NULL; - Mat Gmn=NULL; - NodeSets* nodesets=NULL; - Vec ys=NULL; - Vec ys0=NULL; - - /*intermediary: */ - IoModel* iomodel=NULL; - int sub_analysis_type; - - _printf_(" fill model with matlab workspace data\n"); - iomodel=new IoModel(MODEL); - - _printf_(" specifying analysis\n"); - if (analysis_type!=0){ - iomodel->analysis_type=analysis_type; - } - /*sub_analysis_type should default to none, as it was not specified: */ - sub_analysis_type=NoneAnalysisEnum; - - if (sub_analysis_type!=0){ - iomodel->sub_analysis_type=sub_analysis_type; - } - - _printf_(" create datasets:\n"); - CreateDataSets(&elements,&nodes,&vertices,&materials,&constraints,&loads,¶meters,iomodel,MODEL); - - _printf_(" create degrees of freedom: \n"); - Dofx( &partition,&tpartition,elements,nodes, vertices,parameters); - - _printf_(" create single point constraints: \n"); - SpcNodesx( &yg, nodes,constraints); - - _printf_(" create rigid body constraints:\n"); - MpcNodesx( &Rmg, nodes,constraints); - - _printf_(" create node sets:\n"); - BuildNodeSetsx(&nodesets, nodes); - - _printf_(" reducing single point constraints vector:\n"); - Reducevectorgtosx(&ys,yg->vector,nodesets); - - _printf_(" normalizing rigid body constraints matrix:\n"); - NormalizeConstraintsx(&Gmn, Rmg,nodesets); - - _printf_(" configuring element and loads:\n"); - ConfigureObjectsx(elements, loads, nodes, vertices, materials,parameters); - - _printf_(" process parameters:\n"); - ProcessParamsx( parameters, partition->vector); - - _printf_(" free ressources:\n"); - delete iomodel; - - /*Use all the data created to create an femmodel: */ - femmodel=new FemModel(elements,nodes,vertices, constraints,loads,materials,parameters, - partition,tpartition,yg,Rmg,Gmn,nodesets,ys,ys0); - - /*Add to femmodels in model: */ - femmodels->AddObject((Object*)femmodel); - -} -/*}}}1*/ -/*FUNCTION Model::DeepEcho {{{1*/ - -void Model::DeepEcho(void){ - - femmodels->Echo(); - -} -/*}}}1*/ -/*FUNCTION Model::Echo {{{1*/ - -void Model::Echo(void){ - - printf("Models echo: \n"); - printf(" number of fem models: %i\n",femmodels->Size()); - -} -/*}}}1*/ - -/*FindParam: */ -/*FUNCTION Model::FindParam(bool* pparameter,int enum_type {{{1*/ - -int Model::FindParam(bool* pparameter,int enum_type){ - - FemModel* femmodel=NULL; - - /*no analysis_type or sub_analysis_type, find the parameter in the first dataset we find: */ - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*recover first femmodel: */ - femmodel=(FemModel*)femmodels->GetObjectByOffset(0); - - if(!femmodel)return 0; - - femmodel->parameters->FindParam(pparameter,enum_type); - -} -/*}}}1*/ -/*FUNCTION Model::FindParam(int* pparameter,int enum_type {{{1*/ - -int Model::FindParam(int* pparameter,int enum_type){ - - FemModel* femmodel=NULL; - - /*no analysis_type or sub_analysis_type, find the parameter in the first dataset we find: */ - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*recover first femmodel: */ - femmodel=(FemModel*)femmodels->GetObjectByOffset(0); - - if(!femmodel)return 0; - - femmodel->parameters->FindParam(pparameter,enum_type); - -} -/*}}}1*/ -/*FUNCTION Model::FindParam(double* pparameter,int enum_type ){{{1*/ -int Model::FindParam(double* pparameter,int enum_type){ - - FemModel* femmodel=NULL; - - /*no analysis_type or sub_analysis_type, find the parameter in the first dataset we find: */ - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*recover first femmodel: */ - femmodel=(FemModel*)femmodels->GetObjectByOffset(0); - - if(!femmodel)return 0; - - femmodel->parameters->FindParam(pparameter,enum_type); - - -} -/*}}}1*/ -/*FUNCTION Model::FindParam(double** pparameter,int* pM, int *pN,int enum_type) {{{1*/ -int Model::FindParam(double** pparameter,int* pM, int *pN,int enum_type){ - - FemModel* femmodel=NULL; - - /*no analysis_type or sub_analysis_type, find the parameter in the first dataset we find: */ - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*recover first femmodel: */ - femmodel=(FemModel*)femmodels->GetObjectByOffset(0); - - if(!femmodel)return 0; - - femmodel->parameters->FindParam(pparameter,pM, pN,enum_type); - - -} -/*}}}1*/ -/*FUNCTION Model::FindParam(double** pparameter,int* pM, int enum_type) {{{1*/ -int Model::FindParam(double** pparameter,int* pM, int enum_type){ - - FemModel* femmodel=NULL; - - /*no analysis_type or sub_analysis_type, find the parameter in the first dataset we find: */ - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*recover first femmodel: */ - femmodel=(FemModel*)femmodels->GetObjectByOffset(0); - - if(!femmodel)return 0; - - femmodel->parameters->FindParam(pparameter,pM, enum_type); - - -} -/*}}}1*/ -/*FUNCTION Model::FindParam(char** pparameter,int enum_type) {{{1*/ -int Model::FindParam(char** pparameter,int enum_type){ - - FemModel* femmodel=NULL; - - /*no analysis_type or sub_analysis_type, find the parameter in the first dataset we find: */ - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*recover first femmodel: */ - femmodel=(FemModel*)femmodels->GetObjectByOffset(0); - - if(!femmodel)return 0; - - femmodel->parameters->FindParam(pparameter,enum_type); - -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysisAndSub(bool* pparameter,int enum_type,int analysis_type,int sub_analysis_type) {{{1*/ -int Model::FindParamByAnalysisAndSub(bool* pparameter,int enum_type,int analysis_type,int sub_analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type,sub_analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysisAndSub(int* pparameter,int enum_type,int analysis_type,int sub_analysis_type) {{{1*/ -int Model::FindParamByAnalysisAndSub(int* pparameter,int enum_type,int analysis_type,int sub_analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type,sub_analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysisAndSub(double* pparameter,int enum_type,int analysis_type,int sub_analysis_type) {{{1*/ -int Model::FindParamByAnalysisAndSub(double* pparameter,int enum_type,int analysis_type,int sub_analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type,sub_analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysisAndSub(double** pparameter,int* pM, int* pN,int enum_type,int analysis_type,int sub_analysis_type) {{{1*/ -int Model::FindParamByAnalysisAndSub(double** pparameter,int* pM, int* pN,int enum_type,int analysis_type,int sub_analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type,sub_analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,pM, pN,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysisAndSub(double** pparameter,int* pM, int enum_type,int analysis_type,int sub_analysis_type) {{{1*/ -int Model::FindParamByAnalysisAndSub(double** pparameter,int* pM, int enum_type,int analysis_type,int sub_analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type,sub_analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,pM, enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysisAndSub(char** pparameter,int enum_type,int analysis_type,int sub_analysis_type) {{{1*/ -int Model::FindParamByAnalysisAndSub(char** pparameter,int enum_type,int analysis_type,int sub_analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type,sub_analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysis(bool* pparameter,int enum_type,int analysis_type) {{{1*/ -int Model::FindParamByAnalysis(bool* pparameter,int enum_type,int analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysis(int* pparameter,int enum_type,int analysis_type) {{{1*/ -int Model::FindParamByAnalysis(int* pparameter,int enum_type,int analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysis(double* pparameter,int enum_type,int analysis_type) {{{1*/ -int Model::FindParamByAnalysis(double* pparameter,int enum_type,int analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysis(double** pparameter,int* pM, int* pN,int enum_type,int analysis_type) {{{1*/ -int Model::FindParamByAnalysis(double** pparameter,int* pM, int* pN,int enum_type,int analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,pM,pN,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysis(double** pparameter,int* pM, int enum_type,int analysis_type) {{{1*/ -int Model::FindParamByAnalysis(double** pparameter,int* pM, int enum_type,int analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,pM,enum_type); -} -/*}}}1*/ -/*FUNCTION Model::FindParamByAnalysis(char** pparameter,int enum_type,int analysis_type) {{{1*/ -int Model::FindParamByAnalysis(char** pparameter,int enum_type,int analysis_type){ - - FemModel* femmodel=NULL; - - /*find the correct formulation: */ - femmodel=this->GetFormulation(analysis_type); - - if(!femmodel)return 0; - - /*extract our parameter from the found formulation: */ - femmodel->parameters->FindParam(pparameter,enum_type); -} -/*}}}1*/ - -/*Update:*/ -/*FUNCTION Model::UpdateInputsFromConstant(double constant, int name);{{{1*/ -void Model::UpdateInputsFromConstant(double constant, int name){ - - int i; - FemModel* femmodel=NULL; - - for(i=0;ifemmodels->Size();i++){ - femmodel=(FemModel*)this->femmodels->GetObjectByOffset(i); - UpdateInputsFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,constant,name); - } - -} -/*}}}1*/ -/*FUNCTION Model::UpdateInputsFromConstant(int constant, int name);{{{1*/ -void Model::UpdateInputsFromConstant(int constant, int name){ - - int i; - FemModel* femmodel=NULL; - - for(i=0;ifemmodels->Size();i++){ - femmodel=(FemModel*)this->femmodels->GetObjectByOffset(i); - UpdateInputsFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,constant,name); - } - -} -/*}}}1*/ -/*FUNCTION Model::UpdateInputsFromConstant(bool constant, int name);{{{1*/ -void Model::UpdateInputsFromConstant(bool constant, int name){ - - int i; - FemModel* femmodel=NULL; - - for(i=0;ifemmodels->Size();i++){ - femmodel=(FemModel*)this->femmodels->GetObjectByOffset(i); - UpdateInputsFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,constant,name); - } - -} -/*}}}1*/ -/*FUNCTION Model::UpdateInputsFromVector(double* vector,int name, int type);{{{1*/ -void Model::UpdateInputsFromVector(double* vector,int name, int type){ - - int i; - FemModel* femmodel=NULL; - - if(vector==NULL)return; //don't bother - - for(i=0;ifemmodels->Size();i++){ - femmodel=(FemModel*)this->femmodels->GetObjectByOffset(i); - UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vector,name,type); - } - -} -/*}}}1*/ -/*FUNCTION Model::UpdateInputsFromVector(Vec vector,int name, int type);{{{1*/ -void Model::UpdateInputsFromVector(Vec vector,int name, int type){ - - int i; - FemModel* femmodel=NULL; - double* serial_vector=NULL; - - if(vector==NULL)return; //don't bother - - VecToMPISerial(&serial_vector,vector); - - //call double* idential routine. - this->UpdateInputsFromVector(serial_vector,name,type); - - /*Free ressources:*/ - xfree((void**)&serial_vector); -} -/*}}}1*/ -/*FUNCTION Model::UpdateInputsFromSolution(double* vector,int analysis_type, int sub_analysis_type);{{{1*/ -void Model::UpdateInputsFromSolution(double* vector,int analysis_type, int sub_analysis_type){ - - int i; - FemModel* femmodel=NULL; - - for(i=0;ifemmodels->Size();i++){ - femmodel=(FemModel*)this->femmodels->GetObjectByOffset(i); - UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vector,analysis_type,sub_analysis_type); - } - -} -/*}}}1*/ -/*FUNCTION Model::UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters, double* qmu_part,int qmu_npart);{{{1:*/ -void Model::UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters, double* qmu_part,int qmu_npart){ - ISSMERROR("not supported yet!"); -} -/*}}}*/ - - -/*diverse: */ -/*FUNCTION Model::GetActiveFormulation {{{1*/ -FemModel* Model::GetActiveFormulation(){return active;} -/*}}}1*/ -/*FUNCTION Model::GetFormulation(int analysis_type,int sub_analysis_type) {{{1*/ - -FemModel* Model::GetFormulation(int analysis_type,int sub_analysis_type){ - - int i; - FemModel* femmodel=NULL; - int femmodel_analysis_type; - int femmodel_sub_analysis_type; - int found=0; - - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*find femmodel with correct analysis_type and sub_analysis_type: */ - for(i=0;iSize();i++){ - - femmodel=(FemModel*)femmodels->GetObjectByOffset(i); - femmodel->parameters->FindParam(&femmodel_analysis_type,AnalysisTypeEnum); - femmodel->parameters->FindParam(&femmodel_sub_analysis_type,SubAnalysisTypeEnum); - - if((analysis_type==femmodel_analysis_type) && (sub_analysis_type==femmodel_sub_analysis_type)){ - found=1; - break; - } - } - - if(!found)return NULL; - else return femmodel; -} -/*}}}1*/ -/*FUNCTION Model::GetFormulation(int analysis_type) {{{1*/ - -FemModel* Model::GetFormulation(int analysis_type){ - - int i; - FemModel* femmodel=NULL; - int femmodel_analysis_type; - int femmodel_sub_analysis_type; - int found=0; - - if (!femmodels->Size())ISSMERROR(" no fem models were found!"); - - /*find femmodel with correct analysis_type and sub_analysis_type: */ - for(i=0;iSize();i++){ - - femmodel=(FemModel*)femmodels->GetObjectByOffset(i); - femmodel->parameters->FindParam(&femmodel_analysis_type,AnalysisTypeEnum); - - if((analysis_type==femmodel_analysis_type)){ - found=1; - break; - } - } - - if(!found)return NULL; - else return femmodel; -} -/*}}}1*/ -/*FUNCTION Model::SetActiveFormulation {{{1*/ -void* Model::SetActiveFormulation(FemModel* femmodel){ - active=femmodel; -} -/*}}}1*/ - Index: sm/trunk/src/c/objects/Model.h =================================================================== --- /issm/trunk/src/c/objects/Model.h (revision 3981) +++ (revision ) @@ -1,76 +1,0 @@ -/*! \file Model.h - * \brief: header file for model object - */ - -#ifndef _MODEL_H_ -#define _MODEL_H_ - -/*Headers:*/ -/*{{{1*/ -#include "../include/include.h" -struct FemModel; -class DataSet; -/*}}}*/ - -class Model{ - - private: - - /*femmodels for each formulation: dynamic, can have as many*/ - DataSet* femmodels; - - /*active femmodel: points to a FemModel in the femmodels dataset*/ - FemModel* active; - - public: - - Model(); - ~Model(); - - void Echo(); - void DeepEcho(); - - void AddFormulation(ConstDataHandle MODEL, int analysis_type,int sub_analysis_type); - void AddFormulation(ConstDataHandle MODEL, int analysis_type); - - /*all overloaded forms of the FindParam routine: */ - int FindParam(bool* pparameter,int enum_type); - int FindParam(int* pparameter,int enum_type); - int FindParam(double* pparameter,int enum_type); - int FindParam(double** pparameter,int* pM,int* pN,int enum_type); - int FindParam(double** pparameter,int* pM,int enum_type); - int FindParam(char** pparameter,int enum_type); - - int FindParamByAnalysisAndSub(bool* pparameter,int enum_type,int analysis_type,int sub_analysis_type); - int FindParamByAnalysisAndSub(int* pparameter,int enum_type,int analysis_type,int sub_analysis_type); - int FindParamByAnalysisAndSub(double* pparameter,int enum_type,int analysis_type,int sub_analysis_type); - int FindParamByAnalysisAndSub(double** pparameter,int* pM, int enum_type,int analysis_type,int sub_analysis_type); - int FindParamByAnalysisAndSub(double** pparameter,int* pM, int* pN,int enum_type,int analysis_type,int sub_analysis_type); - int FindParamByAnalysisAndSub(char** pparameter,int enum_type,int analysis_type,int sub_analysis_type); - - int FindParamByAnalysis(bool* pparameter,int enum_type,int analysis_type); - int FindParamByAnalysis(int* pparameter,int enum_type,int analysis_type); - int FindParamByAnalysis(double* pparameter,int enum_type,int analysis_type); - int FindParamByAnalysis(double** pparameter,int* pM,int enum_type,int analysis_type); - int FindParamByAnalysis(double** pparameter,int* pM,int* pN,int enum_type,int analysis_type); - int FindParamByAnalysis(char** pparameter,int enum_type,int analysis_type); - - FemModel* GetFormulation(int analysis_type,int sub_analysis_type); - FemModel* GetFormulation(int analysis_type); - - void UpdateInputsFromVector(Vec vector, int name, int type); - void UpdateInputsFromVector(double* vector, int name, int type); - void UpdateInputsFromVector(int* vector, int name, int type); - void UpdateInputsFromVector(bool* vector, int name, int type); - void UpdateInputsFromConstant(double constant, int name); - void UpdateInputsFromConstant(int constant, int name); - void UpdateInputsFromConstant(bool constant, int name); - void UpdateInputsFromSolution(double* solution, int analysis_type, int sub_analysis_type); - void UpdateFromDakota(double* variables,char** variables_descriptors,int numvariables,DataSet* parameters, double* qmu_part,int qmu_npart); - - - FemModel* GetActiveFormulation(); - void* SetActiveFormulation(FemModel* femmodel); - -}; -#endif /* _MODEL_H */ Index: /issm/trunk/src/c/solutions/diagnostic.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic.cpp (revision 3981) +++ /issm/trunk/src/c/solutions/diagnostic.cpp (revision 3982) @@ -27,6 +27,6 @@ bool control_analysis=false; - /*Model: */ - Model* model=NULL; + /*FemModel: */ + FemModel* femmodel=NULL; /*Results: */ @@ -59,31 +59,31 @@ lockname=argv[4]; - /*Initialize model structure: */ + /*Initialize femmodel structure: */ MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); - model=new Model(); + femmodel=new FemModel(5); /*Open handle to data on disk: */ fid=pfopen(inputfilename,"rb"); - _printf_("read and create finite element model:\n"); - _printf_("\n reading diagnostic horiz model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum); + _printf_("read and create finite element femmodel:\n"); + _printf_("\n reading diagnostic horiz femmodel data:\n"); + femmodel->AddAnalysis(fid,DiagnosticHorizAnalysisEnum); - _printf_("\n reading diagnostic vert model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum); + _printf_("\n reading diagnostic vert femmodel data:\n"); + femmodel->AddAnalysis(fid,DiagnosticVertAnalysisEnum); - _printf_("\n reading diagnostic stokes model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum); + _printf_("\n reading diagnostic stokes femmodel data:\n"); + femmodel->AddAnalysis(fid,DiagnosticStokesAnalysisEnum); - _printf_("\n reading diagnostic hutter model data:\n"); - model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum); + _printf_("\n reading diagnostic hutter femmodel data:\n"); + femmodel->AddAnalysis(fid,DiagnosticHutterAnalysisEnum); - _printf_("\n reading surface and bed slope computation model data:\n"); - model->AddFormulation(fid,SlopecomputeAnalysisEnum); + _printf_("\n reading surface and bed slope computation femmodel data:\n"); + femmodel->AddAnalysis(fid,SlopecomputeAnalysisEnum); /*get parameters: */ - model->FindParam(&qmu_analysis,QmuAnalysisEnum); - model->FindParam(&control_analysis,ControlAnalysisEnum); - model->FindParam(&waitonlock,WaitOnLockEnum); + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum); MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); @@ -95,5 +95,5 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=diagnostic_core(model); + results=diagnostic_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); @@ -103,5 +103,5 @@ _printf_("call computational core:\n"); MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - results=control_core(model); + results=control_core(femmodel); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); @@ -117,5 +117,5 @@ #ifdef _HAVE_DAKOTA_ MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); - Qmux(model,DiagnosticAnalysisEnum,NoneAnalysisEnum); + Qmux(femmodel,DiagnosticAnalysisEnum,NoneAnalysisEnum); MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); #else @@ -130,10 +130,10 @@ /*Free ressources */ - delete model; + delete femmodel; delete results; /*Get finish time and close*/ MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( ); - _printf_("\n %-34s %f seconds \n","Model initialization elapsed time:",finish_init-start_init); + _printf_("\n %-34s %f seconds \n","FemModel initialization elapsed time:",finish_init-start_init); _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core); _printf_("\n %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60); Index: /issm/trunk/src/c/solutions/diagnostic_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 3981) +++ /issm/trunk/src/c/solutions/diagnostic_core.cpp (revision 3982) @@ -11,15 +11,8 @@ #include "../include/include.h" -Results* diagnostic_core(Model* model){ +Results* diagnostic_core(FemModel* femmodel){ extern int my_rank; int dummy; - - /*fem models: */ - FemModel* fem_dh=NULL; - FemModel* fem_dv=NULL; - FemModel* fem_dhu=NULL; - FemModel* fem_ds=NULL; - FemModel* fem_sl=NULL; /*output: */ @@ -70,37 +63,30 @@ results=new Results(); - //first recover parameters common to all solutions - model->FindParam(&verbose,VerboseEnum); - model->FindParam(&dim,DimEnum); - model->FindParam(&ishutter,IsHutterEnum); - model->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); - model->FindParam(&numberofnodes,NumberOfNodesEnum); - model->FindParam(&isstokes,IsStokesEnum); - model->FindParam(&stokesreconditioning,StokesReconditioningEnum); - model->FindParam(&numrifts,NumRiftsEnum); - model->FindParam(&qmu_analysis,QmuAnalysisEnum); - - /*recover fem models: */ - fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum); - fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum); - fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum); - fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum); - fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum); + //first recover parameters needed to drive the solution + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&ishutter,IsHutterEnum); + femmodel->parameters->FindParam(&ismacayealpattyn,IsMacAyealPattynEnum); + femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum); + femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum); + femmodel->parameters->FindParam(&numrifts,NumRiftsEnum); + femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum); //specific parameters for specific models - fem_dh->parameters->FindParam(&numberofdofspernode_dh,NumberOfDofsPerNodeEnum); - fem_sl->parameters->FindParam(&numberofdofspernode_sl,NumberOfDofsPerNodeEnum); - fem_ds->parameters->FindParam(&numberofdofspernode_ds,NumberOfDofsPerNodeEnum); + femmodel->parameters->FindParam(&numberofdofspernode_dh,NumberOfDofsPerNodeEnum); + femmodel->parameters->FindParam(&numberofdofspernode_sl,NumberOfDofsPerNodeEnum); + femmodel->parameters->FindParam(&numberofdofspernode_ds,NumberOfDofsPerNodeEnum); //for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: */ if(qmu_analysis){ - fem_dh->parameters->FindParam(&vx,&dummy,VxEnum); model->UpdateInputsFromVector(vx,VxEnum,VertexEnum); - fem_dh->parameters->FindParam(&vy,&dummy,VyEnum); model->UpdateInputsFromVector(vy,VyEnum,VertexEnum); - fem_dh->parameters->FindParam(&vz,&dummy,VzEnum); model->UpdateInputsFromVector(vz,VzEnum,VertexEnum); + femmodel->parameters->FindParam(&vx,&dummy,VxEnum); femmodel->UpdateInputsFromVector(vx,VxEnum,VertexEnum); + femmodel->parameters->FindParam(&vy,&dummy,VyEnum); femmodel->UpdateInputsFromVector(vy,VyEnum,VertexEnum); + femmodel->parameters->FindParam(&vz,&dummy,VzEnum); femmodel->UpdateInputsFromVector(vz,VzEnum,VertexEnum); } /*Compute slopes: */ - slope_core(&surfaceslopex,&surfaceslopey,fem_sl,SurfaceAnalysisEnum); - slope_core(&bedslopex,&bedslopey,fem_sl,BedAnalysisEnum); + slope_core(&surfaceslopex,&surfaceslopey,femmodel,SurfaceAnalysisEnum); + slope_core(&bedslopex,&bedslopey,femmodel,BedAnalysisEnum); /*Update: */ Index: /issm/trunk/src/c/solutions/slope_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/slope_core.cpp (revision 3981) +++ /issm/trunk/src/c/solutions/slope_core.cpp (revision 3982) @@ -9,5 +9,5 @@ #include "../modules/modules.h" -void slope_core(Vec* pslopex,Vec* pslopey,FemModel* fem, int AnalysisEnum){ +void slope_core(Vec* pslopex,Vec* pslopey,FemModel* femmodel, int AnalysisEnum){ /*parameters: */ @@ -24,8 +24,8 @@ /*Recover some parameters: */ - fem->parameters->FindParam(&verbose,VerboseEnum); - fem->parameters->FindParam(&dim,DimEnum); - fem->parameters->FindParam(&isstokes,IsStokesEnum); - fem->parameters->FindParam(&ishutter,IsHutterEnum); + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&dim,DimEnum); + femmodel->parameters->FindParam(&isstokes,IsStokesEnum); + femmodel->parameters->FindParam(&ishutter,IsHutterEnum); if(verbose)_printf_("%s\n","computing slope (x and y derivatives)..."); @@ -59,12 +59,12 @@ /*Call on core computations: */ - diagnostic_core_linear(&slopex,fem,SlopecomputeAnalysisEnum,xanalysis); - diagnostic_core_linear(&slopey,fem,SlopecomputeAnalysisEnum,yanalysis); + diagnostic_core_linear(&slopex,femmodel,SlopecomputeAnalysisEnum,xanalysis); + diagnostic_core_linear(&slopey,femmodel,SlopecomputeAnalysisEnum,yanalysis); /*extrude if we are in 3D: */ if (dim==3){ if(verbose)_printf_("%s\n","extruding slopes in 3d..."); - FieldExtrudex( slopex, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,"slopex",0); - FieldExtrudex( slopey, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,"slopey",0); + FieldExtrudex( slopex, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"slopex",0); + FieldExtrudex( slopey, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"slopey",0); }