Index: /issm/trunk/src/c/BamgConvertMeshx/BamgConvertMeshx.cpp =================================================================== --- /issm/trunk/src/c/BamgConvertMeshx/BamgConvertMeshx.cpp (revision 3354) +++ /issm/trunk/src/c/BamgConvertMeshx/BamgConvertMeshx.cpp (revision 3354) @@ -0,0 +1,38 @@ +/*!\file BamgConvertMeshx + */ + +#include "./BamgConvertMeshx.h" + +#include "../shared/shared.h" +#include "../include/macros.h" +#include "../toolkits/toolkits.h" + +#include "../Bamgx/objects/BamgObjects.h" + +using namespace bamg; +using namespace std; + +int BamgConvertMeshx(BamgMesh* bamgmesh,BamgGeom* bamggeom,double* index,double* x,double* y,int nods,int nels){ + + /*Intermediary*/ + int i,j,k; + int verbose=0; + int noerr=1; + + /*Options*/ + BamgOpts bamgopts; + + // read mesh + if(verbose) printf("Reading mesh\n"); + Triangles Th(index,x,y,nods,nels); + + //write mesh and geometry + if (verbose) printf("Write Geometry\n"); + Th.Gh.WriteGeometry(bamggeom,&bamgopts); + if (verbose) printf("Write Mesh\n"); + Th.WriteMesh(bamgmesh,&bamgopts); + + /*No error return*/ + return noerr; + +} Index: /issm/trunk/src/c/BamgConvertMeshx/BamgConvertMeshx.h =================================================================== --- /issm/trunk/src/c/BamgConvertMeshx/BamgConvertMeshx.h (revision 3354) +++ /issm/trunk/src/c/BamgConvertMeshx/BamgConvertMeshx.h (revision 3354) @@ -0,0 +1,13 @@ +/*!\file: BamgConvertMeshx.h + * \brief header file for Bamg module + */ + +#ifndef _BAMGCONVERTMESHX_H +#define _BAMGCONVERTMESHX_H + +#include "../objects/objects.h" + +/* local prototypes: */ +int BamgConvertMeshx(BamgMesh* bamgmesh,BamgGeom* bamggeom,double* index,double* x,double* y,int nods,int nels); + +#endif Index: /issm/trunk/src/c/Bamgx/objects/Triangles.cpp =================================================================== --- /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3353) +++ /issm/trunk/src/c/Bamgx/objects/Triangles.cpp (revision 3354) @@ -1307,5 +1307,5 @@ //else st[k]>=0 -> the edge already exist, check else if(st[k]>=0) { - //check that it is not an edge on boundary (should not alrady exist) + //check that it is not an edge on boundary (should not already exist) if (triangles[i].TriangleAdj(j) || triangles[st[k]/3].TriangleAdj((int) (st[k]%3))){ ISSMERROR(exprintf("problem in Geometry reconstruction: an edge on boundary is duplicated (double element?)")); Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 3353) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 3354) @@ -41,6 +41,7 @@ //prognostic int PrognosticAnalysisEnum(void){ return 250; } -int BalancedthicknessAnalysisEnum(void){ return 251; } -int BalancedvelocitiesAnalysisEnum(void){ return 252; } +int Prognostic2AnalysisEnum(void){ return 251; } +int BalancedthicknessAnalysisEnum(void){ return 252; } +int BalancedvelocitiesAnalysisEnum(void){ return 253; } //melting int MeltingAnalysisEnum(void){ return 260; } Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3353) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 3354) @@ -42,4 +42,5 @@ //prognostic int PrognosticAnalysisEnum(void); +int Prognostic2AnalysisEnum(void); int BalancedthicknessAnalysisEnum(void); int BalancedvelocitiesAnalysisEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 3353) +++ /issm/trunk/src/c/Makefile.am (revision 3354) @@ -224,4 +224,8 @@ ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\ + ./ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp\ + ./ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp\ + ./ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp\ + ./ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp\ ./ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp\ ./ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp\ @@ -387,5 +391,7 @@ ./Bamgx/shared/Norm.h \ ./Bamgx/shared/OppositeAngle.h \ - ./Bamgx/shared/shared.h + ./Bamgx/shared/shared.h\ + ./BamgConvertMeshx/BamgConvertMeshx.cpp\ + ./BamgConvertMeshx/BamgConvertMeshx.h libISSM_a_CXXFLAGS = -fPIC -DMATLAB -D_SERIAL_ -ansi -D_GNU_SOURCE -fno-omit-frame-pointer -pthread -D_CPP_ $(CXXOPTFLAGS) @@ -599,4 +605,8 @@ ./ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ ./ModelProcessorx/Prognostic/CreateParametersPrognostic.cpp\ + ./ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp\ + ./ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp\ + ./ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp\ + ./ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp\ ./ModelProcessorx/Balancedthickness/CreateElementsNodesAndMaterialsBalancedthickness.cpp\ ./ModelProcessorx/Balancedthickness/CreateConstraintsBalancedthickness.cpp\ @@ -777,5 +787,7 @@ ./Bamgx/shared/Norm.h \ ./Bamgx/shared/OppositeAngle.h \ - ./Bamgx/shared/shared.h + ./Bamgx/shared/shared.h\ + ./BamgConvertMeshx/BamgConvertMeshx.cpp\ + ./BamgConvertMeshx/BamgConvertMeshx.h libpISSM_a_CXXFLAGS = -fPIC -D_PARALLEL_ -D_C_ $(CXXOPTFLAGS) Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 3353) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 3354) @@ -84,4 +84,12 @@ } + else if (iomodel->analysis_type==Prognostic2AnalysisEnum()){ + + CreateElementsNodesAndMaterialsPrognostic2(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateConstraintsPrognostic2(pconstraints,iomodel,iomodel_handle); + CreateLoadsPrognostic2(ploads,iomodel,iomodel_handle); + CreateParametersPrognostic2(pparameters,iomodel,iomodel_handle); + + } else if (iomodel->analysis_type==BalancedthicknessAnalysisEnum()){ Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3353) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 3354) @@ -248,4 +248,10 @@ void CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + /*prognostic2:*/ + void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateConstraintsPrognostic2(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersPrognostic2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + /*balancedthickness:*/ void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp (revision 3354) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateConstraintsPrognostic2.cpp (revision 3354) @@ -0,0 +1,76 @@ +/* + * CreateConstraintsPrognostic2.c: + */ + + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../IoModel.h" + + +void CreateConstraintsPrognostic2(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + + int i; + int count; + + DataSet* constraints = NULL; + Spc* spc = NULL; + + /*spc intermediary data: */ + int spc_sid; + int spc_node; + int spc_dof; + double spc_value; + + double* spcthickness=NULL; + + /*Create constraints: */ + constraints = new DataSet(ConstraintsEnum()); + + /*Fetch data: */ + IoModelFetchData(&spcthickness,NULL,NULL,iomodel_handle,"spcthickness"); + + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((iomodel->my_grids[i]==1)){ + #endif + + if ((int)spcthickness[2*i]){ + + /*This grid needs to be spc'd: */ + + spc_sid=count; + spc_node=i+1; + spc_dof=1; //we enforce first translation degree of freedom, for temperature + spc_value=*(spcthickness+2*i+1); + + spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); + constraints->AddObject(spc); + count++; + } + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + constraints->Presort(); + + /*Free data: */ + xfree((void**)&spcthickness); + + cleanup_and_return: + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3354) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateElementsNodesAndMaterialsPrognostic2.cpp (revision 3354) @@ -0,0 +1,535 @@ +/* + * CreateElementsNodesAndMaterialsPrognostic2.c: + */ + + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../IoModel.h" + + +void CreateElementsNodesAndMaterialsPrognostic2(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + int i,j,k,n; + extern int my_rank; + extern int num_procs; + + /*DataSets: */ + DataSet* elements = NULL; + DataSet* nodes = NULL; + DataSet* materials = NULL; + + /*Objects: */ + Node* node = NULL; + Tria* tria = NULL; + Penta* penta = NULL; + Matice* matice = NULL; + Matpar* matpar = NULL; + + /*output: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int* my_grids=NULL; + double* my_bordergrids=NULL; + + /*intermediary: */ + int elements_width; //size of elements + double B_avg; + + /*tria constructor input: */ + int tria_id; + int tria_mid; + int tria_mparid; + int tria_numparid; + int tria_g[3]; + double tria_h[3]; + double tria_s[3]; + double tria_b[3]; + double tria_k[3]; + double tria_melting[3]; + double tria_accumulation[3]; + double tria_geothermalflux[3]; + int tria_friction_type; + double tria_p; + double tria_q; + int tria_shelf; + bool tria_onwater; + + /*matice constructor input: */ + int matice_mid; + double matice_B; + double matice_n; + + /*penta constructor input: */ + + int penta_id; + int penta_mid; + int penta_mparid; + int penta_numparid; + int penta_g[6]; + double penta_h[6]; + double penta_s[6]; + double penta_b[6]; + double penta_k[6]; + int penta_friction_type; + double penta_p; + double penta_q; + int penta_shelf; + int penta_onbed; + int penta_onsurface; + int penta_collapse; + double penta_melting[6]; + double penta_accumulation[6]; + double penta_geothermalflux[6]; + int penta_thermal_steadystate; + bool penta_onwater; + + /*matpar constructor input: */ + int matpar_mid; + double matpar_rho_ice; + double matpar_rho_water; + double matpar_heatcapacity; + double matpar_thermalconductivity; + double matpar_latentheat; + double matpar_beta; + double matpar_meltingpoint; + double matpar_mixed_layer_capacity; + double matpar_thermal_exchange_velocity; + double matpar_g; + + /* node constructor input: */ + int node_id; + int node_partitionborder=0; + double node_x[3]; + double node_sigma; + int node_onbed; + int node_onsurface; + int node_onshelf; + int node_onsheet; + int node_upper_node_id; + int node_numdofs; + + + #ifdef _PARALLEL_ + /*Metis partitioning: */ + int range; + Vec gridborder=NULL; + int my_numgrids; + int* all_numgrids=NULL; + int gridcount; + int count; + #endif + int first_grid_index; + + /*First create the elements, nodes and material properties: */ + elements = new DataSet(ElementsEnum()); + nodes = new DataSet(NodesEnum()); + materials = new DataSet(MaterialsEnum()); + + /*Width of elements: */ + if(strcmp(iomodel->meshtype,"2d")==0){ + elements_width=3; //tria elements + } + else{ + ISSMERROR("not implemented yet!"); + 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->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements2d); + + /*Used later on: */ + my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); + #endif + + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + + /*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->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + + for (i=0;inumberofelements;i++){ + + #ifdef _PARALLEL_ + /*!All elements have been partitioned above, only create elements for this CPU: */ + if(my_rank==epart[i]){ + #endif + + + /*ids: */ + tria_id=i+1; //matlab indexing. + tria_mid=-1; //no need for materials + tria_mparid=-1; //no need for materials + tria_numparid=1; + + /*vertices offsets: */ + tria_g[0]=(int)*(iomodel->elements+elements_width*i+0); + tria_g[1]=(int)*(iomodel->elements+elements_width*i+1); + tria_g[2]=(int)*(iomodel->elements+elements_width*i+2); + + /*thickness,surface and bed:*/ + tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. + tria_h[1]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+1)-1)); + tria_h[2]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+2)-1)) ; + + tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1)); + tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1)); + tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1)); + + tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1)); + tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1)); + tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1)); + + /*element on iceshelf?:*/ + tria_shelf=(int)*(iomodel->elementoniceshelf+i); + tria_onwater=(bool)*(iomodel->elementonwater+i); + + /*Create tria element using its constructor:*/ + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_numparid,tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_onwater); + + /*Add tria element to elements dataset: */ + elements->AddObject(tria); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonwater); + + } + else{ // if (strcmp(meshtype,"2d")==0) + ISSMERROR(exprintf("not implemented yet")); + + /*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->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + + for (i=0;inumberofelements;i++){ + #ifdef _PARALLEL_ + /*We are using our element partition to decide which elements will be created on this node: */ + if(my_rank==epart[i]){ + #endif + + + /*name and id: */ + penta_id=i+1; //matlab indexing. + penta_mid=-1; + penta_mparid=-1; //no need for materials + penta_numparid=1; + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(iomodel->elements+elements_width*i+j); + penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + } + + /*diverse: */ + penta_shelf=(int)*(iomodel->elementoniceshelf+i); + penta_onbed=(int)*(iomodel->elementonbed+i); + penta_onsurface=(int)*(iomodel->elementonsurface+i); + penta_collapse=1; + penta_onwater=(bool)*(iomodel->elementonwater+i); + + + /*Create Penta using its constructor:*/ + penta= new Penta( penta_id,penta_mid,penta_mparid,penta_numparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux, + penta_thermal_steadystate,penta_onwater); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elementonwater); + + } //if (strcmp(meshtype,"2d")==0) + + #ifdef _PARALLEL_ + /*From the element partitioning, we can determine which grids are on the inside of this cpu's + *element partition, and which are on its border with other nodes:*/ + gridborder=NewVec(iomodel->numberofnodes); + + for (i=0;inumberofnodes;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + #ifdef _ISSM_DEBUG_ + VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); + #endif + + VecToMPISerial(&my_bordergrids,gridborder); + + #ifdef _ISSM_DEBUG_ + if(my_rank==0){ + for (i=0;inumberofnodes;i++){ + printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); + } + } + #endif + #endif + + /*Add one constant material property to materials: */ + matpar_mid=1; //put it at the end of the materials + matpar_g=iomodel->g; + matpar_rho_ice=iomodel->rho_ice; + matpar_rho_water=iomodel->rho_water; + matpar_thermalconductivity=iomodel->thermalconductivity; + matpar_heatcapacity=iomodel->heatcapacity; + matpar_latentheat=iomodel->latentheat; + matpar_beta=iomodel->beta; + matpar_meltingpoint=iomodel->meltingpoint; + matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; + matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; + + /*Create matpar object using its constructor: */ + matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, + matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, + matpar_thermal_exchange_velocity,matpar_g); + + /*Add to materials datset: */ + materials->AddObject(matpar); + + /*Partition penalties in 3d: */ + if(strcmp(iomodel->meshtype,"3d")==0){ + + /*Get penalties: */ + IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); + + if(iomodel->numpenalties){ + + iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int)); + #ifdef _SERIAL_ + for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=1; + #else + for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=-1; + + for(i=0;inumpenalties;i++){ + first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1); + if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids + /*All grids that are being penalised belong to this node's internal grid partition.:*/ + iomodel->penaltypartitioning[i]=1; + } + if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border + iomodel->penaltypartitioning[i]=0; + } + } + #endif + } + + /*Free penalties: */ + xfree((void**)&iomodel->penalties); + } + + /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. + We can therefore determine which grids are internal to this node's partition + and which ones are shared with other nodes because they are on the border of this node's partition. Knowing + that, go and create the grids*/ + + /*Create nodes from x,y,z, as well as the spc values on those grids: */ + + /*First fetch data: */ + if (strcmp(iomodel->meshtype,"3d")==0){ + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); + } + 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->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); + IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + + + /*Get number of dofs per node: */ + DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((my_grids[i]==1)){ + #endif + + node_id=i+1; //matlab indexing + + + + #ifdef _PARALLEL_ + if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border + node_partitionborder=1; + } + else{ + node_partitionborder=0; + } + #else + node_partitionborder=0; + #endif + + node_x[0]=iomodel->x[i]; + node_x[1]=iomodel->y[i]; + node_x[2]=iomodel->z[i]; + node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); + + node_onbed=(int)iomodel->gridonbed[i]; + node_onsurface=(int)iomodel->gridonsurface[i]; + node_onshelf=(int)iomodel->gridoniceshelf[i]; + node_onsheet=(int)iomodel->gridonicesheet[i]; + + if (strcmp(iomodel->meshtype,"3d")==0){ + if (isnan(iomodel->uppernodes[i])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)iomodel->uppernodes[i]; + } + } + else{ + /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ + node_upper_node_id=node_id; + } + + /*Create node using its constructor: */ + node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + + /*set single point constraints.: */ + if (strcmp(iomodel->meshtype,"3d")==0){ + /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ + if (!iomodel->gridonbed[i]){ + for(k=1;k<=node_numdofs;k++){ + node->FreezeDof(k); + } + } + } + /*Add node to nodes dataset: */ + nodes->AddObject(node); + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + elements->Presort(); + nodes->Presort(); + materials->Presort(); + + /*Clean fetched data: */ + 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->uppernodes); + xfree((void**)&iomodel->gridonicesheet); + xfree((void**)&iomodel->gridoniceshelf); + + + /*Keep partitioning information into iomodel*/ + iomodel->epart=epart; + iomodel->my_grids=my_grids; + iomodel->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + +} Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp (revision 3354) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateLoadsPrognostic2.cpp (revision 3354) @@ -0,0 +1,33 @@ +/*! \file CreateLoadsPrognostic2.c: + */ + + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../IoModel.h" + + +void CreateLoadsPrognostic2(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + DataSet* loads = NULL; + + /*Create loads: */ + loads = new DataSet(LoadsEnum()); + + + /*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: */ + loads->Presort(); + + cleanup_and_return: + + /*Assign output pointer: */ + *ploads=loads; + +} + + Index: /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp (revision 3354) +++ /issm/trunk/src/c/ModelProcessorx/Prognostic2/CreateParametersPrognostic2.cpp (revision 3354) @@ -0,0 +1,155 @@ +/*!\file: CreateParametersPrognostic2.cpp + * \brief driver for creating parameters dataset, for prognostic analysis. + */ + + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../IoModel.h" + +void CreateParametersPrognostic2(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ + + Param* param = NULL; + DataSet* parameters=NULL; + int count; + int i; + int dim; + + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + double* u_g=NULL; + double* pressure=NULL; + double* temperature=NULL; + double* thickness=NULL; + double* surface=NULL; + double* bed=NULL; + double* accumulation=NULL; + double* melting=NULL; + + /*recover parameters : */ + parameters=*pparameters; + + count=parameters->Size(); + + /*Get vx and vy: */ + IoModelFetchData(&vx,NULL,NULL,iomodel_handle,"vx"); + IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy"); + IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz"); + + u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double)); + + if(vx) for(i=0;inumberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts; + if(vy) for(i=0;inumberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts; + if(vz) for(i=0;inumberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts; + + count++; + param= new Param(count,"u_g",DOUBLEVEC); + param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3); + parameters->AddObject(param); + + + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + xfree((void**)&u_g); + + /*Get pressure if 3d iomodel: */ + parameters->FindParam(&dim,"dim"); + if (dim==3){ + IoModelFetchData(&pressure,NULL,NULL,iomodel_handle,"pressure"); + + count++; + param= new Param(count,"p_g",DOUBLEVEC); + if(pressure) param->SetDoubleVec(pressure,iomodel->numberofnodes,1); + else param->SetDoubleVec(pressure,0,0); + parameters->AddObject(param); + + /*Free pressure: */ + xfree((void**)&pressure); + } + + /*Get temperature if 3d iomodel: */ + parameters->FindParam(&dim,"dim"); + if (dim==3){ + IoModelFetchData(&temperature,NULL,NULL,iomodel_handle,"temperature"); + + count++; + param= new Param(count,"t_g",DOUBLEVEC); + if(temperature) param->SetDoubleVec(temperature,iomodel->numberofnodes,1); + else param->SetDoubleVec(temperature,0,0); + parameters->AddObject(param); + + /*Free temperature: */ + xfree((void**)&temperature); + } + + /*Get thickness: */ + IoModelFetchData(&thickness,NULL,NULL,iomodel_handle,"thickness"); + + count++; + param= new Param(count,"h_g",DOUBLEVEC); + if(thickness) param->SetDoubleVec(thickness,iomodel->numberofnodes,1); + else param->SetDoubleVec(thickness,0,0); + parameters->AddObject(param); + + /*Free thickness: */ + xfree((void**)&thickness); + + /*Get surface: */ + IoModelFetchData(&surface,NULL,NULL,iomodel_handle,"surface"); + + count++; + param= new Param(count,"s_g",DOUBLEVEC); + if(surface) param->SetDoubleVec(surface,iomodel->numberofnodes,1); + else param->SetDoubleVec(surface,0,0); + parameters->AddObject(param); + + /*Free surface: */ + xfree((void**)&surface); + + /*Get bed: */ + IoModelFetchData(&bed,NULL,NULL,iomodel_handle,"bed"); + + count++; + param= new Param(count,"b_g",DOUBLEVEC); + if(bed) param->SetDoubleVec(bed,iomodel->numberofnodes,1); + else param->SetDoubleVec(bed,0,0); + parameters->AddObject(param); + + /*Free bed: */ + xfree((void**)&bed); + + /*Get melting: */ + IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting"); + if(melting) for(i=0;inumberofnodes;i++)melting[i]=melting[i]/iomodel->yts; + + count++; + param= new Param(count,"m_g",DOUBLEVEC); + if(melting) param->SetDoubleVec(melting,iomodel->numberofnodes,1); + else param->SetDoubleVec(melting,0,1); + parameters->AddObject(param); + + /*Free melting: */ + xfree((void**)&melting); + + /*Get accumulation: */ + IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation"); + if(accumulation) for(i=0;inumberofnodes;i++)accumulation[i]=accumulation[i]/iomodel->yts; + + count++; + param= new Param(count,"a_g",DOUBLEVEC); + if(accumulation) param->SetDoubleVec(accumulation,iomodel->numberofnodes,1); + else param->SetDoubleVec(accumulation,0,0); + parameters->AddObject(param); + + /*Free accumulation: */ + xfree((void**)&accumulation); + + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk/src/c/issm.h =================================================================== --- /issm/trunk/src/c/issm.h (revision 3353) +++ /issm/trunk/src/c/issm.h (revision 3354) @@ -65,4 +65,5 @@ #include "./MassFluxx/MassFluxx.h" #include "./Bamgx/Bamgx.h" +#include "./BamgConvertMeshx/BamgConvertMeshx.h" Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 3353) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 3354) @@ -53,4 +53,7 @@ numdofs=1; } + else if (analysis_type==Prognostic2AnalysisEnum()){ + numdofs=1; + } else if (analysis_type==BalancedthicknessAnalysisEnum()){ numdofs=1; Index: /issm/trunk/src/m/classes/public/bamg.m =================================================================== --- /issm/trunk/src/m/classes/public/bamg.m (revision 3353) +++ /issm/trunk/src/m/classes/public/bamg.m (revision 3354) @@ -285,5 +285,5 @@ %Fill in rest of fields: md.type='2d'; -md.numberofelements=length(md.elements); +md.numberofelements=size(md.elements,1); md.numberofgrids=length(md.x); md.z=zeros(md.numberofgrids,1); Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m =================================================================== --- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 3353) +++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m (revision 3354) @@ -106,5 +106,5 @@ %SIZE NUMBEROFELEMENTS fields={'elements'}; -checklength(md,fields,md.numberofelements); +checksize(md,fields,[md.numberofelements 3]); %SIZE NUMBEROFGRIDS Index: /issm/trunk/src/m/classes/public/mesh/meshconvert.m =================================================================== --- /issm/trunk/src/m/classes/public/mesh/meshconvert.m (revision 3354) +++ /issm/trunk/src/m/classes/public/mesh/meshconvert.m (revision 3354) @@ -0,0 +1,47 @@ +function md=meshconvert(md,varargin) +%CONVERTMESH - convert mesh to bamg mesh +% +% Usage: +% md=meshconvert(md); +% md=meshconvert(md,index,x,y); + +if nargin~=1 & nargin~=4, + help meshconvert + error('meshconvert error message: bad usage'); +end + +if nargin==1, + x=md.x; + y=md.y; + index=md.elements; +else + x=varargin{1}; + y=varargin{2}; + index=varargin{3}; +end + +%call Bamg +[bamgmesh_out bamggeom_out]=BamgConvertMesh(index,x,y); + +% plug results onto model +md.bamg=struct(); +md.bamg.mesh=bamgmesh(bamgmesh_out); +md.bamg.geometry=bamggeom(bamggeom_out); +md.x=bamgmesh_out.Vertices(:,1); +md.y=bamgmesh_out.Vertices(:,2); +md.elements=bamgmesh_out.Triangles(:,1:3); +md.segments=bamgmesh_out.Segments(:,1:3); +md.segmentmarkers=bamgmesh_out.Segments(:,4); + +%Fill in rest of fields: +md.type='2d'; +md.numberofelements=size(md.elements,1); +md.numberofgrids=length(md.x); +md.z=zeros(md.numberofgrids,1); +md.gridonbed=ones(md.numberofgrids,1); +md.gridonwater=zeros(md.numberofgrids,1); +md.gridonsurface=ones(md.numberofgrids,1); +md.elementonbed=ones(md.numberofelements,1); +md.elementonsurface=ones(md.numberofelements,1); +md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1; +md.counter=1; Index: /issm/trunk/src/m/classes/public/process_solve_options.m =================================================================== --- /issm/trunk/src/m/classes/public/process_solve_options.m (revision 3353) +++ /issm/trunk/src/m/classes/public/process_solve_options.m (revision 3354) @@ -20,5 +20,5 @@ %check solution type is supported -if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient',... +if ~ismemberi(analysis_type,{'diagnostic','prognostic','prognostic2','thermal','steadystate','parameters','transient',... 'balancedthickness','balancedvelocities','slopecompute'}), error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']); Index: /issm/trunk/src/m/classes/public/solve.m =================================================================== --- /issm/trunk/src/m/classes/public/solve.m (revision 3353) +++ /issm/trunk/src/m/classes/public/solve.m (revision 3354) @@ -60,5 +60,8 @@ elseif md.analysis_type==PrognosticAnalysisEnum, - md=prognostic(md);; + md=prognostic(md); + +elseif md.analysis_type==Prognostic2AnalysisEnum, + md=prognostic2(md); elseif md.analysis_type==ThermalAnalysisEnum, Index: /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m =================================================================== --- /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m (revision 3353) +++ /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m (revision 3354) @@ -85,4 +85,8 @@ end +if enum==Prognostic2AnalysisEnum(), + string='prognostic2'; +end + if enum==BalancedthicknessAnalysisEnum(), string='balancedthickness'; Index: /issm/trunk/src/m/enum/BalancedthicknessAnalysisEnum.m =================================================================== --- /issm/trunk/src/m/enum/BalancedthicknessAnalysisEnum.m (revision 3353) +++ /issm/trunk/src/m/enum/BalancedthicknessAnalysisEnum.m (revision 3354) @@ -7,3 +7,3 @@ % macro=BalancedthicknessAnalysisEnum() -macro=251; +macro=252; Index: /issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m =================================================================== --- /issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m (revision 3353) +++ /issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m (revision 3354) @@ -7,3 +7,3 @@ % macro=BalancedvelocitiesAnalysisEnum() -macro=252; +macro=253; Index: /issm/trunk/src/m/enum/Prognostic2AnalysisEnum.m =================================================================== --- /issm/trunk/src/m/enum/Prognostic2AnalysisEnum.m (revision 3354) +++ /issm/trunk/src/m/enum/Prognostic2AnalysisEnum.m (revision 3354) @@ -0,0 +1,9 @@ +function macro=Prognostic2AnalysisEnum() +%PROGNOSTIC2ANALYSISENUM - Enum of Prognostic2Analysis +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=Prognostic2AnalysisEnum() + +macro=251; Index: sm/trunk/src/m/enum/SteadyAnalysisEnum.m =================================================================== --- /issm/trunk/src/m/enum/SteadyAnalysisEnum.m (revision 3353) +++ (revision ) @@ -1,9 +1,0 @@ -function macro=TransientAnalysisEnum() -%TRANSIENTANALYSISENUM - Enum of TransientAnalysis -% -% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum -% -% Usage: -% macro=TransientAnalysisEnum() - -macro=231; Index: /issm/trunk/src/m/solutions/jpl/prognostic2.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2.m (revision 3354) +++ /issm/trunk/src/m/solutions/jpl/prognostic2.m (revision 3354) @@ -0,0 +1,38 @@ +function md=prognostic2(md) +%PROGNOSITC2 - prognostic2 solution sequence. +% +% Usage: +% md=prognostic2(md) + + %timing + t1=clock; + + %Build all models requested for diagnostic simulation + models.analysis_type=Prognostic2AnalysisEnum; %needed for processresults + + displaystring(md.verbose,'%s',['reading prognostic2 model data']); + md.analysis_type=Prognostic2AnalysisEnum; models.p=CreateFemModel(md); + error('STOP here'); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + displaystring(md.verbose,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes); + inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double'); + + displaystring(md.verbose,'\n%s',['call computational core:']); + results=prognostic2_core(models,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum()); + + displaystring(md.verbose,'\n%s',['load results...']); + if ~isstruct(md.results), md.results=struct(); end + md.results.prognostic2=processresults(models, results); + + %stop timing + t2=clock; + displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/jpl/prognostic2_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 3354) +++ /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 3354) @@ -0,0 +1,24 @@ +function results=prognostic2_core(models,inputs,analysis_type,sub_analysis_type) +%PROGNOSTIC2_CORE - linear solution sequence +% +% Usage: +% h_g=prognostic2_core(m,inputs,analysis_type,sub_analysis_type) + + %get FE model + m=models.p; + results.time=0; + results.step=1; + + displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']); + %Take only the first two dofs of m.parameters.u_g + u_g=get(inputs,'velocity',[1 1 0 0]); + u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity'); + inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes); + + displaystring(m.parameters.verbose,'\n%s',['call computational core:']); + results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type); + + displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']); + results.h_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0); + +end %end function Index: /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.cpp =================================================================== --- /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.cpp (revision 3354) +++ /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.cpp (revision 3354) @@ -0,0 +1,82 @@ +/*\file BamgConvertMesh.c + *\brief: bamg module. + */ +#include "./BamgConvertMesh.h" + +void mexFunction( int nlhs, mxArray* plhs[], int nrhs, const mxArray* prhs[]){ + + /*input: */ + double* index=NULL; + int index_rows; + double* x=NULL; + int x_cols; + double* y=NULL; + int y_rows; + int y_cols; + + /*Output*/ + BamgMesh bamgmesh; + BamgGeom bamggeom; + mxArray* bamgmesh_mat=NULL; + mxArray* bamggeom_mat=NULL; + + /*Intermediary*/ + int nods; + int nels; + int verbose=2; + + /*Boot module: */ + MODULEBOOT(); + + /*checks on arguments on the matlab side: */ + CheckNumMatlabArguments(nlhs,NLHS,nrhs,NRHS,__FUNCT__,&BamgConvertMeshUsage); + + /*Input datasets: */ + if (verbose) printf("Fetching inputs\n"); + FetchData(&index,&nels,&index_rows,INDEXHANDLE); + FetchData(&x,&nods,&x_cols,XHANDLE); + FetchData(&y,&y_rows,&y_cols,YHANDLE); + + /*Check inputs*/ + if (nels<0){ + ISSMERROR("Number of elements must be positive, check index number of lines"); + } + if (nods<0){ + ISSMERROR("Number of nods must be positive, check x and y sizes"); + } + if (index_rows!=3){ + ISSMERROR("index should have 3 columns"); + } + if (y_rows!=nods){ + ISSMERROR("x and y do not have the same length"); + } + if (x_cols>1 || y_cols>1){ + ISSMERROR("x and y should have only one column"); + } + + /* Run core computations: */ + if (verbose) printf("Call core\n"); + BamgConvertMeshx(&bamgmesh,&bamggeom,index,x,y,nods,nels); + + /*Generate output Matlab Structures*/ + WriteData(&bamgmesh_mat,&bamgmesh); + WriteData(&bamggeom_mat,&bamggeom); + + /*assign output datasets: */ + *BAMGMESHOUT=bamgmesh_mat; + *BAMGGEOMOUT=bamggeom_mat; + + /*end module: */ + MODULEEND(); +} + +void BamgConvertMeshUsage(void) +{ + _printf_("BAMGCONVERTMESH - convert [x y index] to a bamg geom and mesh geom"); + _printf_("\n"); + _printf_(" Usage:\n"); + _printf_(" [bamggeom bamgmesh]=BamgConvertMesh(index,x,y);\n"); + _printf_(" index: index of the mesh\n"); + _printf_(" x,y: coordinates of the nodes\n"); + _printf_("\n"); +} Index: /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.h =================================================================== --- /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.h (revision 3354) +++ /issm/trunk/src/mex/BamgConvertMesh/BamgConvertMesh.h (revision 3354) @@ -0,0 +1,32 @@ +/*!\file BamgConvertMesh.h + * \brief: prototype for Data Interpolation mex module. + */ + +#ifndef _BAMGCONVERTMESH_H +#define _BAMGCONVERTMESH_H + +/* local prototypes: */ +void BamgConvertMeshUsage(void); + +#include "../../c/issm.h" + +#undef __FUNCT__ +#define __FUNCT__ "BamgConvertMesh" + + +/* serial input macros: */ +#define INDEXHANDLE prhs[0] +#define XHANDLE prhs[1] +#define YHANDLE prhs[2] + +/* serial output macros: */ +#define BAMGMESHOUT (mxArray**)&plhs[0] +#define BAMGGEOMOUT (mxArray**)&plhs[1] + +/* serial arg counts: */ +#undef NLHS +#define NLHS 2 +#undef NRHS +#define NRHS 3 + +#endif Index: /issm/trunk/src/mex/Makefile.am =================================================================== --- /issm/trunk/src/mex/Makefile.am (revision 3353) +++ /issm/trunk/src/mex/Makefile.am (revision 3354) @@ -7,4 +7,5 @@ bin_PROGRAMS = BuildNodeSets\ Bamg\ + BamgConvertMesh\ ComputePressure\ ConfigureObjects \ @@ -88,4 +89,7 @@ Bamg/Bamg.h +BamgConvertMesh_SOURCES = BamgConvertMesh/BamgConvertMesh.cpp\ + BamgConvertMesh/BamgConvertMesh.h + ComputePressure_SOURCES = ComputePressure/ComputePressure.cpp\ ComputePressure/ComputePressure.h