Index: /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/Control/CreateParametersControl.cpp (revision 117) @@ -0,0 +1,131 @@ +/*!\file: CreateParametersControl.cpp + * \brief driver for creating parameters dataset, for control analysis. + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateParameters" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + +void CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount){ + + int i; + + Param* param = NULL; + int count; + int analysis_type; + int numberofdofspernode; + + double* fit=NULL; + double* optscal=NULL; + double* maxiter=NULL; + double* control_parameter=NULL; + + /*Get counter : */ + count=*pcount; + + /*control_type: */ + count++; + param= new Param(count,"control_type",STRING); + param->SetString(model->control_type); + parameters->AddObject(param); + + + /*nsteps: */ + count++; + param= new Param(count,"nsteps",INTEGER); + param->SetInteger(model->nsteps); + parameters->AddObject(param); + + /*tolx: */ + count++; + param= new Param(count,"tolx",DOUBLE); + param->SetDouble(model->tolx); + parameters->AddObject(param); + + /*mincontrolconstraint: */ + count++; + param= new Param(count,"mincontrolconstraint",DOUBLE); + param->SetDouble(model->mincontrolconstraint); + parameters->AddObject(param); + + /*maxcontrolconstraint: */ + count++; + param= new Param(count,"maxcontrolconstraint",DOUBLE); + param->SetDouble(model->maxcontrolconstraint); + parameters->AddObject(param); + + /*epsvel: */ + count++; + param= new Param(count,"epsvel",DOUBLE); + param->SetDouble(model->epsvel); + parameters->AddObject(param); + + /*meanvel: */ + count++; + param= new Param(count,"meanvel",DOUBLE); + param->SetDouble(model->meanvel); + parameters->AddObject(param); + + /*Now, recover fit, optscal and maxiter as vectors: */ + ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat"); + ModelFetchData((void**)&model->optscal,NULL,NULL,model_handle,"optscal","Matrix","Mat"); + ModelFetchData((void**)&model->maxiter,NULL,NULL,model_handle,"maxiter","Matrix","Mat"); + + count++; + param= new Param(count,"fit",DOUBLEVEC); + param->SetDoubleVec(model->fit,model->nsteps); + parameters->AddObject(param); + + count++; + param= new Param(count,"optscal",DOUBLEVEC); + param->SetDoubleVec(model->optscal,model->nsteps); + parameters->AddObject(param); + + count++; + param= new Param(count,"maxiter",DOUBLEVEC); + param->SetDoubleVec(model->maxiter,model->nsteps); + parameters->AddObject(param); + + xfree((void**)&model->fit); + xfree((void**)&model->optscal); + xfree((void**)&model->maxiter); + + /*Get vx_obs, vy_obs, and the parameter value: */ + ModelFetchData((void**)&model->vx_obs,NULL,NULL,model_handle,"vx_obs","Matrix","Mat"); + ModelFetchData((void**)&model->vy_obs,NULL,NULL,model_handle,"vy_obs","Matrix","Mat"); + ModelFetchData((void**)&control_parameter,NULL,NULL,model_handle,model->control_type,"Matrix","Mat"); + + for(i=0;inumberofnodes;i++)model->vx_obs[i]=model->vx_obs[i]/model->yts; + for(i=0;inumberofnodes;i++)model->vy_obs[i]=model->vy_obs[i]/model->yts; + + count++; + param= new Param(count,"vx_obs",DOUBLEVEC); + param->SetDoubleVec(model->vx_obs,model->numberofnodes); + parameters->AddObject(param); + + count++; + param= new Param(count,"vy_obs",DOUBLEVEC); + param->SetDoubleVec(model->vy_obs,model->numberofnodes); + parameters->AddObject(param); + + count++; + param= new Param(count,"control_parameter",DOUBLEVEC); + param->SetDoubleVec(control_parameter,model->numberofnodes); + parameters->AddObject(param); + + xfree((void**)&model->vx_obs); + xfree((void**)&model->vy_obs); + xfree((void**)&control_parameter); + + + /*Assign output pointer: */ + *pcount=count; +} + + Index: /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/CreateConstraints.cpp (revision 117) @@ -0,0 +1,41 @@ +/*!\file: CreateConstraints.cpp + * \brief general driver for creating constraints + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "CreateConstraints" + +#include "./Model.h" +#include "../shared/shared.h" + +void CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ + + /*This is just a high level driver: */ + if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ + CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle); + } + /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){ + CreateConstraintsDiagnosticBaseVert(pconstraints,model,model_handle); + } + else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ + CreateConstraintsDiagnosticVert(pconstraints,model,model_handle); + } + else if (strcmp(model->analysis_type,"melting")==0){ + CreateConstraintsMelting(pconstraints,model,model_handle); + } + else if (strcmp(model->analysis_type,"prognostic")==0){ + CreateConstraintsPrognostic(pconstraints,model,model_handle); + } + else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ + CreateConstraintsThermal(pconstraints,model,model_handle); + }*/ + else{ + throw ErrorException(__FUNCT__," analysis_type not supported yet!"); + } +} Index: sm/trunk/src/c/ModelProcessorx/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateConstraintsDiagnosticHoriz.cpp (revision 116) +++ (revision ) @@ -1,90 +1,0 @@ -/* - * CreateConstraintsDiagnosticHoriz.c: - */ - -#undef __FUNCT__ -#define __FUNCT__ "CreateConstraintsDiagnosticHoriz" - -#include "../DataSet/DataSet.h" -#include "../toolkits/toolkits.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "./ModelProcessorx.h" - - -int CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ - - - int noerr=1; - 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* dirichletvalues_diag=NULL; - double* gridondirichlet_diag=NULL; - - /*Create constraints: */ - constraints = new DataSet(ConstraintsEnum()); - - /*Fetch data: */ - ModelFetchData((void**)&gridondirichlet_diag,NULL,NULL,model_handle,"gridondirichlet_diag","Matrix","Mat"); - ModelFetchData((void**)&dirichletvalues_diag,NULL,NULL,model_handle,"dirichletvalues_diag","Matrix","Mat"); - - 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((model->my_grids[i]==1)){ - #endif - - if ((int)gridondirichlet_diag[i]){ - - /*This grid needs to be spc'd to vx_obs and vy_obs:*/ - - spc_sid=count; - spc_node=i+1; - spc_dof=1; //we enforce first x translation degree of freedom - spc_value=*(dirichletvalues_diag+2*i+0)/model->yts; - - spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); - constraints->AddObject(spc); - count++; - - spc_sid=count; - spc_node=i+1; - spc_dof=2; //we enforce first y translation degree of freedom - spc_value=*(dirichletvalues_diag+2*i+1)/model->yts; - - 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**)&gridondirichlet_diag); - xfree((void**)&dirichletvalues_diag); - - /*Assign output pointer: */ - *pconstraints=constraints; - return noerr; -} Index: /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterials.cpp (revision 117) @@ -0,0 +1,50 @@ +/*!\file: CreateElementsNodesAndMaterials.cpp + * \brief general driver for creating elements, nodes and materials + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterials" + +#include "./Model.h" +#include "../shared/shared.h" + + +void CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ + + /*This is just a high level driver: */ + if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ + + CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle); + + } + /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){ + + CreateElementsNodesAndMaterialsDataSetsDiagnosticBaseVert(pelements,pnodes,pmaterials, model,model->analysis_type); + + } + else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ + + CreateElementsNodesAndMaterialsDataSetsDiagnosticVert(pelements,pnodes,pmaterials, model,model->analysis_type); + } + else if (strcmp(model->analysis_type,"melting")==0){ + + CreateElementsNodesAndMaterialsDataSetsMelting(pelements,pnodes,pmaterials, model,model->analysis_type); + } + else if (strcmp(model->analysis_type,"prognostic")==0){ + + CreateElementsNodesAndMaterialsDataSetsPrognostic(pelements,pnodes,pmaterials, model,model->analysis_type); + } + else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ + + CreateElementsNodesAndMaterialsDataSetsThermal(pelements,pnodes,pmaterials, model,model->analysis_type); + }*/ + else{ + throw ErrorException(__FUNCT__," analysis_type not supported yet!"); + } +} Index: sm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 116) +++ (revision ) @@ -1,651 +1,0 @@ -/* - * CreateElementsNodesAndMaterialsDiagnosticHoriz.c: - */ - -#undef __FUNCT__ -#define __FUNCT__ "CreateElementsNodesAndMaterialsDiagnosticHoriz" - -#include "../DataSet/DataSet.h" -#include "../toolkits/toolkits.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "./ModelProcessorx.h" -#include "./MeshPartitionx/MeshPartitionx.h" - - -int CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ - - - /*output: int* epart, int* my_grids, double* my_bordergrids*/ - - - int noerr=1; - 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; - - int analysis_type; - - /*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_g[3]; - double tria_h[3]; - double tria_s[3]; - double tria_b[3]; - double tria_k[3]; - int tria_friction_type; - double tria_p; - double tria_q; - int tria_shelf; - double tria_meanvel;/*!scaling ratio for velocities*/ - double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - double tria_viscosity_overshoot; - - /*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_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; - double penta_meanvel;/*!scaling ratio for velocities*/ - double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ - int penta_collapse; - double penta_melting[6]; - double penta_accumulation[6]; - double penta_geothermalflux[6]; - int penta_artdiff; - int penta_thermal_steadystate; - double penta_viscosity_overshoot; - - /*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]; - int node_onbed; - int node_onsurface; - int node_upper_node_id; - int node_numdofs; - - - #ifdef _PARALLEL_ - /*Metis partitioning: */ - int range; - Vec gridborder; - int my_numgrids; - int* all_numgrids=NULL; - int gridcount; - int count; - #endif - int first_grid_index; - - /*Rifts:*/ - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - double* riftpenaltypairs=NULL; - int el1,el2; - - /*Penalty partitioning: */ - int num_grids3d_collapsed; - double* double_penalties_grids3d_collapsed=NULL; - double* double_penalties_grids3d_noncollapsed=NULL; - int grid_ids[6]; - int num_grid_ids; - int grid_id; - - /*Get analysis_type: */ - analysis_type=AnalysisTypeAsEnum(model->analysis_type); - - /*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(model->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(model->meshtype,"2d")==0){ - /*load elements: */ - ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); - } - else{ - /*load elements2d: */ - ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); - } - - MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs); - - /*Free elements and elements2d: */ - xfree((void**)&model->elements); - xfree((void**)&model->elements2d); - - - /*Deal with rifts, they have to be included into one partition only, not several: */ - FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing - el2=(int)*(riftpenaltypairs+7*j+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**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); - - /*Used later on: */ - my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int)); - #endif - - - - /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ - - /*2d mesh: */ - if (strcmp(model->meshtype,"2d")==0){ - - /*Fetch data needed: */ - ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); - ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); - ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); - ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); - ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat"); - ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat"); - ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); - ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); - ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); - ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); - - 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=i+1; //refers to the corresponding material property card - tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card - - /*vertices offsets: */ - tria_g[0]=(int)*(model->elements+elements_width*i+0); - tria_g[1]=(int)*(model->elements+elements_width*i+1); - tria_g[2]=(int)*(model->elements+elements_width*i+2); - - /*thickness,surface and bed:*/ - tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. - tria_h[1]=*(model->thickness+ ((int)*(model->elements+elements_width*i+1)-1)); - tria_h[2]=*(model->thickness+ ((int)*(model->elements+elements_width*i+2)-1)) ; - - tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1)); - tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1)); - tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1)); - - tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1)); - tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1)); - tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1)); - - /*basal drag:*/ - tria_friction_type=(int)model->drag_type; - - tria_k[0]=*(model->drag+ ((int)*(model->elements+elements_width*i+0)-1)); - tria_k[1]=*(model->drag+ ((int)*(model->elements+elements_width*i+1)-1)); - tria_k[2]=*(model->drag+ ((int)*(model->elements+elements_width*i+2)-1)); - - tria_p=model->p[i]; - tria_q=model->q[i]; - - /*element on iceshelf?:*/ - tria_shelf=(int)*(model->elementoniceshelf+i); - - tria_meanvel=model->meanvel; - tria_epsvel=model->epsvel; - - /*viscosity_overshoot*/ - tria_viscosity_overshoot=model->viscosity_overshoot; - - /*Create tria element using its constructor:*/ - tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); - - /*Add tria element to elements dataset: */ - elements->AddObject(tria); - - /*Deal with material property card: */ - matice_mid=i+1; //same as the material id from the geom2 elements. - - /*Average B over 3 grid elements: */ - B_avg=0; - for(j=0;j<3;j++){ - B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1)); - } - B_avg=B_avg/3; - matice_B=B_avg; - matice_n=(double)*(model->n+i); - - /*Create matice using its constructor:*/ - matice= new Matice(matice_mid,matice_B,matice_n); - - /*Add matice element to materials dataset: */ - materials->AddObject(matice); - - #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)*(model->elements+elements_width*i+0)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; - #endif - - #ifdef _PARALLEL_ - }//if(my_rank==epart[i]) - #endif - - }//for (i=0;ielements); - xfree((void**)&model->thickness); - xfree((void**)&model->surface); - xfree((void**)&model->bed); - xfree((void**)&model->drag); - xfree((void**)&model->p); - xfree((void**)&model->q); - xfree((void**)&model->elementoniceshelf); - xfree((void**)&model->B); - xfree((void**)&model->n);*/ - } - else{ // if (strcmp(meshtype,"2d")==0) - - /*Fetch data needed: */ - ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); - ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); - ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); - ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); - ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat"); - ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat"); - ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); - ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); - ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat"); - ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat"); - ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat"); - ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); - ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); - - 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=i+1; //refers to the corresponding material property card - penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card - - /*vertices,thickness,surface,bed and drag: */ - for(j=0;j<6;j++){ - penta_g[j]=(int)*(model->elements+elements_width*i+j); - penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1)); - penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1)); - penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1)); - penta_k[j]=*(model->drag+ ((int)*(model->elements+elements_width*i+j)-1)); - } - - /*basal drag:*/ - penta_friction_type=(int)model->drag_type; - - penta_p=model->p[i]; - penta_q=model->q[i]; - - /*diverse: */ - penta_shelf=(int)*(model->elementoniceshelf+i); - penta_onbed=(int)*(model->elementonbed+i); - penta_onsurface=(int)*(model->elementonsurface+i); - penta_meanvel=model->meanvel; - penta_epsvel=model->epsvel; - - /*viscosity_overshoot*/ - penta_viscosity_overshoot=model->viscosity_overshoot; - - if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base. - penta_collapse=1; - } - else{ - penta_collapse=0; - } - - - /*Create Penta using its constructor:*/ - penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, - penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel, - penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, - penta_thermal_steadystate,penta_viscosity_overshoot); - - /*Add penta element to elements dataset: */ - elements->AddObject(penta); - - - /*Deal with material:*/ - matice_mid=i+1; //same as the material id from the geom2 elements. - /*Average B over 6 element grids: */ - B_avg=0; - for(j=0;j<6;j++){ - B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1)); - } - B_avg=B_avg/6; - matice_B= B_avg; - matice_n=(double)*(model->n+i); - - /*Create matice using its constructor:*/ - matice= new Matice(matice_mid,matice_B,matice_n); - - /*Add matice element to materials dataset: */ - materials->AddObject(matice); - - #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)*(model->elements+elements_width*i+0)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+3)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+4)-1]=1; - my_grids[(int)*(model->elements+elements_width*i+5)-1]=1; - #endif - - #ifdef _PARALLEL_ - }//if(my_rank==epart[i]) - #endif - - }//for (i=0;ielements); - xfree((void**)&model->thickness); - xfree((void**)&model->surface); - xfree((void**)&model->bed); - xfree((void**)&model->drag); - xfree((void**)&model->p); - xfree((void**)&model->q); - xfree((void**)&model->elementoniceshelf); - xfree((void**)&model->elementonbed); - xfree((void**)&model->elementonsurface); - xfree((void**)&model->elements_type); - xfree((void**)&model->n); - xfree((void**)&model->B); - - } //if (strcmp(meshtype,"2d")==0) - - /*Add one constant material property to materials: */ - matpar_mid=model->numberofelements+1; //put it at the end of the materials - matpar_g=model->g; - matpar_rho_ice=model->rho_ice; - matpar_rho_water=model->rho_water; - matpar_thermalconductivity=model->thermalconductivity; - matpar_heatcapacity=model->heatcapacity; - matpar_latentheat=model->latentheat; - matpar_beta=model->beta; - matpar_meltingpoint=model->meltingpoint; - matpar_mixed_layer_capacity=model->mixed_layer_capacity; - matpar_thermal_exchange_velocity=model->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); - - #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(model->numberofnodes); - - for (i=0;inumberofnodes;i++){ - if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); - } - VecAssemblyBegin(gridborder); - VecAssemblyEnd(gridborder); - - #ifdef _DEBUG_ - VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); - #endif - - VecToMPISerial(&my_bordergrids,gridborder); - - #ifdef _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 - - /*Partition penalties in 3d: */ - if(strcmp(model->meshtype,"3d")==0){ - - /*Get penalties: */ - ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); - - if(model->numpenalties){ - - model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int)); - for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=-1; - - for(i=0;inumpenalties;i++){ - first_grid_index=(int)(*(model->penalties+i*model->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.:*/ - model->penaltypartitioning[i]=1; - } - if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border - model->penaltypartitioning[i]=0; - } - } - } - - /*Free penalties: */ - xfree((void**)&model->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(model->meshtype,"3d")==0){ - ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat"); - ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat"); - } - ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat"); - ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat"); - ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat"); - ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat"); - ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); - - - /*Get number of dofs per node: */ - DistributeNumDofs(&node_numdofs,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]=model->x[i]; - node_x[1]=model->y[i]; - node_x[2]=model->z[i]; - - - node_onbed=(int)model->gridonbed[i]; - node_onsurface=(int)model->gridonsurface[i]; - if (strcmp(model->meshtype,"3d")==0){ - if (isnan(model->uppernodes[i])){ - node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. - } - else{ - node_upper_node_id=(int)model->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_onbed,node_onsurface,node_upper_node_id); - - /*set single point constraints.: */ - if (strcmp(model->meshtype,"3d")==0){ - /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ - if (model->deadgrids[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**)&model->deadgrids); - xfree((void**)&model->x); - xfree((void**)&model->y); - xfree((void**)&model->z); - xfree((void**)&model->gridonbed); - xfree((void**)&model->gridonsurface); - xfree((void**)&model->uppernodes); - - cleanup_and_return: - - /*Assign output pointer: */ - *pelements=elements; - *pnodes=nodes; - *pmaterials=materials; - - /*Keep partitioning information into model*/ - model->epart=epart; - model->my_grids=my_grids; - model->my_bordergrids=my_bordergrids; - - /*Free ressources:*/ - #ifdef _PARALLEL_ - xfree((void**)&all_numgrids); - VecFree(&gridborder); - #endif - - return noerr; - -} Index: /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/CreateLoads.cpp (revision 117) @@ -0,0 +1,48 @@ +/*!\file: CreateLoads.cpp + * \brief create loads datasets. General driver. + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#undef __FUNCT__ +#define __FUNCT__ "CreateLoads" + +#include "./Model.h" +#include "../shared/shared.h" + +void CreateLoads(DataSet** ploads, Model* model,ConstDataHandle model_handle){ + + /*This is just a high level driver: */ + if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ + + CreateLoadsDiagnosticHoriz(ploads,model,model_handle); + } + /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){ + + CreateLoadsDiagnosticBaseVert(ploads,model,model_handle); + } + else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ + + CreateLoadsDiagnosticVert(ploads,model,model_handle); + } + else if (strcmp(model->analysis_type,"melting")==0){ + + CreateLoadsMelting(ploads,model,model_handle); + } + else if (strcmp(model->analysis_type,"prognostic")==0){ + + CreateLoadsPrognostic(ploads,model,model_handle); + } + else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ + + CreateLoadsThermal(ploads,model,model_handle); + + }*/ + else{ + throw ErrorException(__FUNCT__," analysis_type not supported yet!"); + } +} Index: sm/trunk/src/c/ModelProcessorx/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateLoadsDiagnosticHoriz.cpp (revision 116) +++ (revision ) @@ -1,281 +1,0 @@ -/*! \file CreateLoadsDiagnosticHoriz.c: - */ - -#undef __FUNCT__ -#define __FUNCT__ "CreateLoadsDiagnosticHoriz" - -#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 "./ModelProcessorx.h" - - -int CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){ - - int noerr=1; - int i,j,counter; - int element; - - extern int my_rank; - extern int num_procs; - - DataSet* loads = NULL; - Icefront* icefront = NULL; - Penpair* penpair = NULL; - - int segment_width; - int element_type; - int i1,i2,i3,i4; - int i0,range; - int grid1,grid2; - - /*rift penpair: */ - double normal[2]; - double length; - double friction; - int fill; - - /*icefront intermediary data: */ - char icefront_type[ICEFRONTSTRING]; - int icefront_element_type; - int icefront_sid; - int icefront_eid; - int icefront_mparid; - int icefront_node_ids[MAX_ICEFRONT_GRIDS]; - double icefront_h[MAX_ICEFRONT_GRIDS]; - double icefront_b[MAX_ICEFRONT_GRIDS]; - - /*penpair intermediary data: */ - int penpair_id; - int penpair_node_ids[2]; - double penpair_penalty_offset; - int penpair_numdofs; - int penpair_dof; - int penpair_penalty_lock; - int penpair_element_ids[2]; - double penpair_friction; - int penpair_fill; - double penpair_normal[2]; - double penpair_length; - - /*Rifts:*/ - int* riftsnumpenaltypairs=NULL; - double** riftspenaltypairs=NULL; - int* riftsfill=NULL; - int* riftsfriction=NULL; - double* riftpenaltypairs=NULL; - int el1,el2; - - /*Create loads: */ - loads = new DataSet(LoadsEnum()); - - /*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids - * referenced by a certain load must belong to the cluster node): */ - ModelFetchData((void**)&model->segmentonneumann_diag,&model->numberofsegs_diag,NULL,model_handle,"segmentonneumann_diag","Matrix","Mat"); - ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat"); - ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); - ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); - - /*First load data:*/ - for (i=0;inumberofsegs_diag;i++){ - - if (strcmp(model->meshtype,"2d")==0){ - segment_width=3; - element_type=TriaEnum(); - } - else{ - segment_width=5; - element_type=PentaEnum(); - } - - - element=(int)(*(model->segmentonneumann_diag+segment_width*i+segment_width-1)-1); //element is in the last column - - #ifdef _PARALLEL_ - if (model->epart[element]!=my_rank){ - /*This load does not belong to this cluster node, as it references an element - *that does not belong to this node's partition. Drop this 'i':*/ - continue; - } - #endif - - icefront_mparid=model->numberofelements+1; //matlab indexing - icefront_sid=i+1; //matlab indexing - icefront_eid=(int)*(model->segmentonneumann_diag+segment_width*i+segment_width-1); //matlab indexing - icefront_element_type=element_type; - - i1=(int)*(model->segmentonneumann_diag+segment_width*i+0); - i2=(int)*(model->segmentonneumann_diag+segment_width*i+1); - - icefront_node_ids[0]=i1; - icefront_node_ids[1]=i2; - - icefront_h[0]=model->thickness[i1-1]; - icefront_h[1]=model->thickness[i2-1]; - - icefront_b[0]=model->bed[i1-1]; - icefront_b[1]=model->bed[i2-1]; - - if ((int)*(model->elements_type+2*element+0)==MacAyealEnum()){ //this is a collapsed 3d element, icefront will be 2d - strcpy(icefront_type,"segment"); - } - else if ((int)*(model->elements_type+2*element+0)==PattynEnum()){ //this is a real 3d element, icefront will be 3d. - strcpy(icefront_type,"quad"); - i3=(int)*(model->segmentonneumann_diag+segment_width*i+2); - i4=(int)*(model->segmentonneumann_diag+segment_width*i+3); - icefront_node_ids[2]=i3; - icefront_node_ids[3]=i4; - - icefront_h[2]=model->thickness[i3-1]; - icefront_h[3]=model->thickness[i4-1]; - - icefront_b[2]=model->bed[i3-1]; - icefront_b[3]=model->bed[i4-1]; - } - else{ - _printf_("%s%s\n",__FUNCT__," error message: element type not supported yet"); - return 0; - } - - icefront = new Icefront(icefront_type,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b); - - loads->AddObject(icefront); - - } - /*Free data: */ - xfree((void**)&model->segmentonneumann_diag); - xfree((void**)&model->elements_type); - xfree((void**)&model->thickness); - xfree((void**)&model->bed); - - - /*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */ - - /*First fetch penalties: */ - if (strcmp(model->meshtype,"3d")==0){ - ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); - } - - counter=0; - if(model->numpenalties){ - - /*First deal with internal grids: */ - if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU - - for(j=0;jnumlayers-1;j++){ - - /*We are pairing grids along a vertical profile.*/ - grid1=(int)*(model->penalties+model->numlayers*i+j); - grid2=(int)*(model->penalties+model->numlayers*i+j+1); - - penpair_id=counter+1; //matlab indexing - penpair_dof=1; - penpair_node_ids[0]=grid1; - penpair_node_ids[1]=grid2; - penpair_penalty_offset=model->penalty_offset; - penpair_numdofs=1; - - penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, - penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); - loads->AddObject(penpair); - - - counter++; - - penpair_id=counter+1; //matlab indexing - penpair_dof=2; - penpair_node_ids[0]=grid1; - penpair_node_ids[1]=grid2; - penpair_penalty_offset=model->penalty_offset; - penpair_numdofs=1; - - penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, - penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); - loads->AddObject(penpair); - - - counter++; - - } - } //for ( i=0; ipenalties); - - /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */ - /*First fetch rifts: */ - if(model->numrifts){ - - for(i=0;inumrifts;i++){ - riftpenaltypairs=model->riftspenaltypairs[i]; - for(j=0;jriftsnumpenaltypairs[i];j++){ - - el1=(int)*(riftpenaltypairs+7*j+2); - #ifdef _PARALLEL_ - if (model->epart[el1-1]!=my_rank){ - /*This load does not belong to this cluster node, as it references an element - *that does not belong to this node's partition. Drop this 'j':*/ - continue; - } - #endif - - /*Ok, retrieve all the data needed to add a penalty between the two grids: */ - el2=(int)*(riftpenaltypairs+7*j+3); - grid1=(int)*(riftpenaltypairs+7*j+0); - grid2=(int)*(riftpenaltypairs+7*j+1); - normal[0]=*(riftpenaltypairs+7*j+4); - normal[1]=*(riftpenaltypairs+7*j+5); - length=*(riftpenaltypairs+7*j+6); - friction=model->riftsfriction[i]; - fill=model->riftsfill[i]; - - penpair_id=counter+1; //matlab indexing - penpair_node_ids[0]=grid1; - penpair_node_ids[1]=grid2; - penpair_element_ids[0]=el1; - penpair_element_ids[1]=el2; - penpair_penalty_offset=model->penalty_offset; - penpair_penalty_lock=model->penalty_lock; - penpair_numdofs=2; - penpair_normal[0]=normal[0]; - penpair_normal[1]=normal[1]; - penpair_length=length; - penpair_friction=friction; - penpair_fill=fill; - - penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, - penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); - - loads->AddObject(penpair); - - counter++; - } - } - } - - /*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(); - - /*Free ressources:*/ - xfree((void**)&riftsnumpenaltypairs); - for(i=0;inumrifts;i++){ - double* temp=riftspenaltypairs[i]; - xfree((void**)&temp); - } - xfree((void**)&riftspenaltypairs); - xfree((void**)&riftsfill); - xfree((void**)&riftsfriction); - - - /*Assign output pointer: */ - *ploads=loads; - - return noerr; -} - - Index: /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/CreateParameters.cpp (revision 117) @@ -0,0 +1,183 @@ +/*!\file: CreateParameters.cpp + * \brief general driver for creating parameters dataset + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateParameters" + +#include "../DataSet/DataSet.h" +#include "../toolkits/toolkits.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "./Model.h" + +void CreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ + + int i; + + DataSet* parameters = NULL; + Param* param = NULL; + int count=1; + int analysis_type; + int numberofdofspernode; + + double* fit=NULL; + double* optscal=NULL; + double* maxiter=NULL; + double* control_parameter=NULL; + + + + /*Initialize dataset: */ + parameters = new DataSet(ParametersEnum()); + + //Get analysis_type: + analysis_type=AnalysisTypeAsEnum(model->analysis_type); + + count++; + param= new Param(count,"analysis_type",INTEGER); + param->SetInteger(analysis_type); + parameters->AddObject(param); + + + /*debug: */ + param= new Param(count,"debug",INTEGER); + param->SetInteger(model->debug); + parameters->AddObject(param); + + /*eps_rel: */ + count++; + param= new Param(count,"eps_rel",DOUBLE); + param->SetDouble(model->eps_rel); + parameters->AddObject(param); + + /*eps_abs: */ + count++; + param= new Param(count,"eps_abs",DOUBLE); + param->SetDouble(model->eps_abs); + parameters->AddObject(param); + + /*yts: */ + count++; + param= new Param(count,"yts",DOUBLE); + param->SetDouble(model->yts); + parameters->AddObject(param); + + /*dt: */ + count++; + param= new Param(count,"dt",DOUBLE); + param->SetDouble(model->dt); + parameters->AddObject(param); + + /*ndt: */ + count++; + param= new Param(count,"ndt",DOUBLE); + param->SetDouble(model->ndt); + parameters->AddObject(param); + + /*penalty_offset: */ + count++; + param= new Param(count,"penalty_offset",DOUBLE); + param->SetDouble(model->penalty_offset); + parameters->AddObject(param); + + /*sparsity: */ + count++; + param= new Param(count,"sparsity",DOUBLE); + param->SetDouble(model->sparsity); + parameters->AddObject(param); + + /*lowmem: */ + count++; + param= new Param(count,"lowmem",INTEGER); + param->SetInteger(model->lowmem); + parameters->AddObject(param); + + /*connectivity: */ + count++; + param= new Param(count,"connectivity",INTEGER); + param->SetInteger(model->connectivity); + parameters->AddObject(param); + + /*beta: */ + count++; + param= new Param(count,"beta",DOUBLE); + param->SetDouble(model->beta); + parameters->AddObject(param); + + /*meltingpoint: */ + count++; + param= new Param(count,"meltingpoint",DOUBLE); + param->SetDouble(model->meltingpoint); + parameters->AddObject(param); + + /*latentheat: */ + count++; + param= new Param(count,"latentheat",DOUBLE); + param->SetDouble(model->latentheat); + parameters->AddObject(param); + + /*heatcapacity: */ + count++; + param= new Param(count,"heatcapacity",DOUBLE); + param->SetDouble(model->heatcapacity); + parameters->AddObject(param); + + /*penalty_melting: */ + count++; + param= new Param(count,"penalty_melting",DOUBLE); + param->SetDouble(model->penalty_melting); + parameters->AddObject(param); + + /*min_thermal_constraints: */ + count++; + param= new Param(count,"min_thermal_constraints",INTEGER); + param->SetInteger(model->min_thermal_constraints); + parameters->AddObject(param); + + /*waitonlock: */ + count++; + param= new Param(count,"waitonlock",INTEGER); + param->SetInteger(model->waitonlock); + parameters->AddObject(param); + + /*solverstring: */ + count++; + param= new Param(count,"solverstring",STRING); + param->SetString(model->solverstring); + parameters->AddObject(param); + + /*plot: */ + count++; + param= new Param(count,"plot",INTEGER); + param->SetInteger(model->plot); + parameters->AddObject(param); + + /*numberofgrids: */ + count++; + param= new Param(count,"numberofnodes",INTEGER); + param->SetInteger(model->numberofnodes); + parameters->AddObject(param); + + /*Deal with numberofdofspernode: */ + DistributeNumDofs(&numberofdofspernode,analysis_type); + + count++; + param= new Param(count,"numberofdofspernode",INTEGER); + param->SetInteger(numberofdofspernode); + parameters->AddObject(param); + + + /*Deal with control parameters: */ + if (analysis_type==ControlAnalysisEnum()){ + CreateParametersControl(parameters,model,model_handle,&count); + } + + /*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + parameters->Presort(); + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 117) @@ -0,0 +1,88 @@ +/* + * CreateConstraintsDiagnosticHoriz.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateConstraintsDiagnosticHoriz" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../Model.h" + + +void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints, Model* model,ConstDataHandle model_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* dirichletvalues_diag=NULL; + double* gridondirichlet_diag=NULL; + + /*Create constraints: */ + constraints = new DataSet(ConstraintsEnum()); + + /*Fetch data: */ + ModelFetchData((void**)&gridondirichlet_diag,NULL,NULL,model_handle,"gridondirichlet_diag","Matrix","Mat"); + ModelFetchData((void**)&dirichletvalues_diag,NULL,NULL,model_handle,"dirichletvalues_diag","Matrix","Mat"); + + 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((model->my_grids[i]==1)){ + #endif + + if ((int)gridondirichlet_diag[i]){ + + /*This grid needs to be spc'd to vx_obs and vy_obs:*/ + + spc_sid=count; + spc_node=i+1; + spc_dof=1; //we enforce first x translation degree of freedom + spc_value=*(dirichletvalues_diag+2*i+0)/model->yts; + + spc = new Spc(spc_sid,spc_node,spc_dof,spc_value); + constraints->AddObject(spc); + count++; + + spc_sid=count; + spc_node=i+1; + spc_dof=2; //we enforce first y translation degree of freedom + spc_value=*(dirichletvalues_diag+2*i+1)/model->yts; + + 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**)&gridondirichlet_diag); + xfree((void**)&dirichletvalues_diag); + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 117) @@ -0,0 +1,648 @@ +/* + * CreateElementsNodesAndMaterialsDiagnosticHoriz.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterialsDiagnosticHoriz" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../Model.h" + + +void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle){ + + + /*output: int* epart, int* my_grids, double* my_bordergrids*/ + + + 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; + + int analysis_type; + + /*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_g[3]; + double tria_h[3]; + double tria_s[3]; + double tria_b[3]; + double tria_k[3]; + int tria_friction_type; + double tria_p; + double tria_q; + int tria_shelf; + double tria_meanvel;/*!scaling ratio for velocities*/ + double tria_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + double tria_viscosity_overshoot; + + /*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_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; + double penta_meanvel;/*!scaling ratio for velocities*/ + double penta_epsvel; /*!minimum velocity to avoid infinite velocity ratios*/ + int penta_collapse; + double penta_melting[6]; + double penta_accumulation[6]; + double penta_geothermalflux[6]; + int penta_artdiff; + int penta_thermal_steadystate; + double penta_viscosity_overshoot; + + /*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]; + int node_onbed; + int node_onsurface; + int node_upper_node_id; + int node_numdofs; + + + #ifdef _PARALLEL_ + /*Metis partitioning: */ + int range; + Vec gridborder; + int my_numgrids; + int* all_numgrids=NULL; + int gridcount; + int count; + #endif + int first_grid_index; + + /*Rifts:*/ + int* riftsnumpenaltypairs; + double** riftspenaltypairs; + int* riftsfill; + double* riftsfriction; + double* riftpenaltypairs=NULL; + int el1,el2; + + /*Penalty partitioning: */ + int num_grids3d_collapsed; + double* double_penalties_grids3d_collapsed=NULL; + double* double_penalties_grids3d_noncollapsed=NULL; + int grid_ids[6]; + int num_grid_ids; + int grid_id; + + /*Get analysis_type: */ + analysis_type=AnalysisTypeAsEnum(model->analysis_type); + + /*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(model->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(model->meshtype,"2d")==0){ + /*load elements: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + } + else{ + /*load elements2d: */ + ModelFetchData((void**)&model->elements2d,NULL,NULL,model_handle,"elements2d","Matrix","Mat"); + } + + MeshPartitionx(&epart, &npart,model->numberofelements,model->numberofnodes,model->elements, model->numberofelements2d,model->numberofnodes2d,model->elements2d,model->numlayers,elements_width, model->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&model->elements); + xfree((void**)&model->elements2d); + + + /*Deal with rifts, they have to be included into one partition only, not several: */ + FetchRifts(&riftsnumpenaltypairs,&riftspenaltypairs,&riftsfill,&riftsfriction,model_handle,model->numrifts); + + for(i=0;inumrifts;i++){ + riftpenaltypairs=model->riftspenaltypairs[i]; + for(j=0;jriftsnumpenaltypairs[i];j++){ + el1=(int)*(riftpenaltypairs+7*j+2)-1; //matlab indexing to c indexing + el2=(int)*(riftpenaltypairs+7*j+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**)&riftsnumpenaltypairs); + for(i=0;inumrifts;i++){ + double* temp=riftspenaltypairs[i]; + xfree((void**)&temp); + } + xfree((void**)&riftspenaltypairs); + xfree((void**)&riftsfill); + xfree((void**)&riftsfriction); + + /*Used later on: */ + my_grids=(int*)xcalloc(model->numberofnodes,sizeof(int)); + #endif + + + + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + + /*2d mesh: */ + if (strcmp(model->meshtype,"2d")==0){ + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat"); + ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat"); + ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); + ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); + ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); + ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); + + 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=i+1; //refers to the corresponding material property card + tria_mparid=model->numberofelements+1;//refers to the corresponding parmat property card + + /*vertices offsets: */ + tria_g[0]=(int)*(model->elements+elements_width*i+0); + tria_g[1]=(int)*(model->elements+elements_width*i+1); + tria_g[2]=(int)*(model->elements+elements_width*i+2); + + /*thickness,surface and bed:*/ + tria_h[0]= *(model->thickness+ ((int)*(model->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. + tria_h[1]=*(model->thickness+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_h[2]=*(model->thickness+ ((int)*(model->elements+elements_width*i+2)-1)) ; + + tria_s[0]=*(model->surface+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_s[1]=*(model->surface+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_s[2]=*(model->surface+ ((int)*(model->elements+elements_width*i+2)-1)); + + tria_b[0]=*(model->bed+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_b[1]=*(model->bed+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_b[2]=*(model->bed+ ((int)*(model->elements+elements_width*i+2)-1)); + + /*basal drag:*/ + tria_friction_type=(int)model->drag_type; + + tria_k[0]=*(model->drag+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_k[1]=*(model->drag+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_k[2]=*(model->drag+ ((int)*(model->elements+elements_width*i+2)-1)); + + tria_p=model->p[i]; + tria_q=model->q[i]; + + /*element on iceshelf?:*/ + tria_shelf=(int)*(model->elementoniceshelf+i); + + tria_meanvel=model->meanvel; + tria_epsvel=model->epsvel; + + /*viscosity_overshoot*/ + tria_viscosity_overshoot=model->viscosity_overshoot; + + /*Create tria element using its constructor:*/ + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_g, tria_h, tria_s, tria_b, tria_k, tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot); + + /*Add tria element to elements dataset: */ + elements->AddObject(tria); + + /*Deal with material property card: */ + matice_mid=i+1; //same as the material id from the geom2 elements. + + /*Average B over 3 grid elements: */ + B_avg=0; + for(j=0;j<3;j++){ + B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1)); + } + B_avg=B_avg/3; + matice_B=B_avg; + matice_n=(double)*(model->n+i); + + /*Create matice using its constructor:*/ + matice= new Matice(matice_mid,matice_B,matice_n); + + /*Add matice element to materials dataset: */ + materials->AddObject(matice); + + #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)*(model->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->thickness); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + xfree((void**)&model->drag); + xfree((void**)&model->p); + xfree((void**)&model->q); + xfree((void**)&model->elementoniceshelf); + xfree((void**)&model->B); + xfree((void**)&model->n);*/ + } + else{ // if (strcmp(meshtype,"2d")==0) + + /*Fetch data needed: */ + ModelFetchData((void**)&model->elements,NULL,NULL,model_handle,"elements","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->surface,NULL,NULL,model_handle,"surface","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + ModelFetchData((void**)&model->drag,NULL,NULL,model_handle,"drag","Matrix","Mat"); + ModelFetchData((void**)&model->p,NULL,NULL,model_handle,"p","Matrix","Mat"); + ModelFetchData((void**)&model->q,NULL,NULL,model_handle,"q","Matrix","Mat"); + ModelFetchData((void**)&model->elementoniceshelf,NULL,NULL,model_handle,"elementoniceshelf","Matrix","Mat"); + ModelFetchData((void**)&model->elementonbed,NULL,NULL,model_handle,"elementonbed","Matrix","Mat"); + ModelFetchData((void**)&model->elementonsurface,NULL,NULL,model_handle,"elementonsurface","Matrix","Mat"); + ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat"); + ModelFetchData((void**)&model->B,NULL,NULL,model_handle,"B","Matrix","Mat"); + ModelFetchData((void**)&model->n,NULL,NULL,model_handle,"n","Matrix","Mat"); + + 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=i+1; //refers to the corresponding material property card + penta_mparid=model->numberofelements+1;//refers to the corresponding parmat property card + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(model->elements+elements_width*i+j); + penta_h[j]=*(model->thickness+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_s[j]=*(model->surface+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_b[j]=*(model->bed+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_k[j]=*(model->drag+ ((int)*(model->elements+elements_width*i+j)-1)); + } + + /*basal drag:*/ + penta_friction_type=(int)model->drag_type; + + penta_p=model->p[i]; + penta_q=model->q[i]; + + /*diverse: */ + penta_shelf=(int)*(model->elementoniceshelf+i); + penta_onbed=(int)*(model->elementonbed+i); + penta_onsurface=(int)*(model->elementonsurface+i); + penta_meanvel=model->meanvel; + penta_epsvel=model->epsvel; + + /*viscosity_overshoot*/ + penta_viscosity_overshoot=model->viscosity_overshoot; + + if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 3 are MacAyeal type Penta. We collapse the formulation on their base. + penta_collapse=1; + } + else{ + penta_collapse=0; + } + + + /*Create Penta using its constructor:*/ + penta= new Penta( penta_id,penta_mid,penta_mparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface,penta_meanvel,penta_epsvel, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux,penta_artdiff, + penta_thermal_steadystate,penta_viscosity_overshoot); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + + /*Deal with material:*/ + matice_mid=i+1; //same as the material id from the geom2 elements. + /*Average B over 6 element grids: */ + B_avg=0; + for(j=0;j<6;j++){ + B_avg+=*(model->B+((int)*(model->elements+elements_width*i+j)-1)); + } + B_avg=B_avg/6; + matice_B= B_avg; + matice_n=(double)*(model->n+i); + + /*Create matice using its constructor:*/ + matice= new Matice(matice_mid,matice_B,matice_n); + + /*Add matice element to materials dataset: */ + materials->AddObject(matice); + + #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)*(model->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+2)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+3)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+4)-1]=1; + my_grids[(int)*(model->elements+elements_width*i+5)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&model->thickness); + xfree((void**)&model->surface); + xfree((void**)&model->bed); + xfree((void**)&model->drag); + xfree((void**)&model->p); + xfree((void**)&model->q); + xfree((void**)&model->elementoniceshelf); + xfree((void**)&model->elementonbed); + xfree((void**)&model->elementonsurface); + xfree((void**)&model->elements_type); + xfree((void**)&model->n); + xfree((void**)&model->B); + + } //if (strcmp(meshtype,"2d")==0) + + /*Add one constant material property to materials: */ + matpar_mid=model->numberofelements+1; //put it at the end of the materials + matpar_g=model->g; + matpar_rho_ice=model->rho_ice; + matpar_rho_water=model->rho_water; + matpar_thermalconductivity=model->thermalconductivity; + matpar_heatcapacity=model->heatcapacity; + matpar_latentheat=model->latentheat; + matpar_beta=model->beta; + matpar_meltingpoint=model->meltingpoint; + matpar_mixed_layer_capacity=model->mixed_layer_capacity; + matpar_thermal_exchange_velocity=model->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); + + #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(model->numberofnodes); + + for (i=0;inumberofnodes;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + #ifdef _DEBUG_ + VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); + #endif + + VecToMPISerial(&my_bordergrids,gridborder); + + #ifdef _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 + + /*Partition penalties in 3d: */ + if(strcmp(model->meshtype,"3d")==0){ + + /*Get penalties: */ + ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); + + if(model->numpenalties){ + + model->penaltypartitioning=(int*)xmalloc(model->numpenalties*sizeof(int)); + for(i=0;inumpenalties;i++)model->penaltypartitioning[i]=-1; + + for(i=0;inumpenalties;i++){ + first_grid_index=(int)(*(model->penalties+i*model->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.:*/ + model->penaltypartitioning[i]=1; + } + if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border + model->penaltypartitioning[i]=0; + } + } + } + + /*Free penalties: */ + xfree((void**)&model->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(model->meshtype,"3d")==0){ + ModelFetchData((void**)&model->deadgrids,NULL,NULL,model_handle,"deadgrids","Matrix","Mat"); + ModelFetchData((void**)&model->uppernodes,NULL,NULL,model_handle,"uppergrids","Matrix","Mat"); + } + ModelFetchData((void**)&model->x,NULL,NULL,model_handle,"x","Matrix","Mat"); + ModelFetchData((void**)&model->y,NULL,NULL,model_handle,"y","Matrix","Mat"); + ModelFetchData((void**)&model->z,NULL,NULL,model_handle,"z","Matrix","Mat"); + ModelFetchData((void**)&model->gridonbed,NULL,NULL,model_handle,"gridonbed","Matrix","Mat"); + ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); + + + /*Get number of dofs per node: */ + DistributeNumDofs(&node_numdofs,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]=model->x[i]; + node_x[1]=model->y[i]; + node_x[2]=model->z[i]; + + + node_onbed=(int)model->gridonbed[i]; + node_onsurface=(int)model->gridonsurface[i]; + if (strcmp(model->meshtype,"3d")==0){ + if (isnan(model->uppernodes[i])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)model->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_onbed,node_onsurface,node_upper_node_id); + + /*set single point constraints.: */ + if (strcmp(model->meshtype,"3d")==0){ + /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ + if (model->deadgrids[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**)&model->deadgrids); + xfree((void**)&model->x); + xfree((void**)&model->y); + xfree((void**)&model->z); + xfree((void**)&model->gridonbed); + xfree((void**)&model->gridonsurface); + xfree((void**)&model->uppernodes); + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + + /*Keep partitioning information into model*/ + model->epart=epart; + model->my_grids=my_grids; + model->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + VecFree(&gridborder); + #endif + +} Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 117) @@ -0,0 +1,278 @@ +/*! \file CreateLoadsDiagnosticHoriz.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateLoadsDiagnosticHoriz" + +#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 "../Model.h" + + +void CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model,ConstDataHandle model_handle){ + + int i,j,counter; + int element; + + extern int my_rank; + extern int num_procs; + + DataSet* loads = NULL; + Icefront* icefront = NULL; + Penpair* penpair = NULL; + + int segment_width; + int element_type; + int i1,i2,i3,i4; + int i0,range; + int grid1,grid2; + + /*rift penpair: */ + double normal[2]; + double length; + double friction; + int fill; + + /*icefront intermediary data: */ + char icefront_type[ICEFRONTSTRING]; + int icefront_element_type; + int icefront_sid; + int icefront_eid; + int icefront_mparid; + int icefront_node_ids[MAX_ICEFRONT_GRIDS]; + double icefront_h[MAX_ICEFRONT_GRIDS]; + double icefront_b[MAX_ICEFRONT_GRIDS]; + + /*penpair intermediary data: */ + int penpair_id; + int penpair_node_ids[2]; + double penpair_penalty_offset; + int penpair_numdofs; + int penpair_dof; + int penpair_penalty_lock; + int penpair_element_ids[2]; + double penpair_friction; + int penpair_fill; + double penpair_normal[2]; + double penpair_length; + + /*Rifts:*/ + int* riftsnumpenaltypairs=NULL; + double** riftspenaltypairs=NULL; + int* riftsfill=NULL; + int* riftsfriction=NULL; + double* riftpenaltypairs=NULL; + int el1,el2; + + /*Create loads: */ + loads = new DataSet(LoadsEnum()); + + /*Create pressure loads as boundary conditions. Pay attention to the partitioning if we are running in parallel (the grids + * referenced by a certain load must belong to the cluster node): */ + ModelFetchData((void**)&model->segmentonneumann_diag,&model->numberofsegs_diag,NULL,model_handle,"segmentonneumann_diag","Matrix","Mat"); + ModelFetchData((void**)&model->elements_type,NULL,NULL,model_handle,"elements_type","Matrix","Mat"); + ModelFetchData((void**)&model->thickness,NULL,NULL,model_handle,"thickness","Matrix","Mat"); + ModelFetchData((void**)&model->bed,NULL,NULL,model_handle,"bed","Matrix","Mat"); + + /*First load data:*/ + for (i=0;inumberofsegs_diag;i++){ + + if (strcmp(model->meshtype,"2d")==0){ + segment_width=3; + element_type=TriaEnum(); + } + else{ + segment_width=5; + element_type=PentaEnum(); + } + + + element=(int)(*(model->segmentonneumann_diag+segment_width*i+segment_width-1)-1); //element is in the last column + + #ifdef _PARALLEL_ + if (model->epart[element]!=my_rank){ + /*This load does not belong to this cluster node, as it references an element + *that does not belong to this node's partition. Drop this 'i':*/ + continue; + } + #endif + + icefront_mparid=model->numberofelements+1; //matlab indexing + icefront_sid=i+1; //matlab indexing + icefront_eid=(int)*(model->segmentonneumann_diag+segment_width*i+segment_width-1); //matlab indexing + icefront_element_type=element_type; + + i1=(int)*(model->segmentonneumann_diag+segment_width*i+0); + i2=(int)*(model->segmentonneumann_diag+segment_width*i+1); + + icefront_node_ids[0]=i1; + icefront_node_ids[1]=i2; + + icefront_h[0]=model->thickness[i1-1]; + icefront_h[1]=model->thickness[i2-1]; + + icefront_b[0]=model->bed[i1-1]; + icefront_b[1]=model->bed[i2-1]; + + if ((int)*(model->elements_type+2*element+0)==MacAyealEnum()){ //this is a collapsed 3d element, icefront will be 2d + strcpy(icefront_type,"segment"); + } + else if ((int)*(model->elements_type+2*element+0)==PattynEnum()){ //this is a real 3d element, icefront will be 3d. + strcpy(icefront_type,"quad"); + i3=(int)*(model->segmentonneumann_diag+segment_width*i+2); + i4=(int)*(model->segmentonneumann_diag+segment_width*i+3); + icefront_node_ids[2]=i3; + icefront_node_ids[3]=i4; + + icefront_h[2]=model->thickness[i3-1]; + icefront_h[3]=model->thickness[i4-1]; + + icefront_b[2]=model->bed[i3-1]; + icefront_b[3]=model->bed[i4-1]; + } + else{ + throw ErrorException(__FUNCT__," element type not supported yet"); + } + + icefront = new Icefront(icefront_type,icefront_sid,icefront_mparid,icefront_eid,icefront_element_type,icefront_node_ids,icefront_h,icefront_b); + + loads->AddObject(icefront); + + } + /*Free data: */ + xfree((void**)&model->segmentonneumann_diag); + xfree((void**)&model->elements_type); + xfree((void**)&model->thickness); + xfree((void**)&model->bed); + + + /*Create penalties loads for linking of collapsed formulation grids with non collapsed grids: */ + + /*First fetch penalties: */ + if (strcmp(model->meshtype,"3d")==0){ + ModelFetchData((void**)&model->penalties,&model->numpenalties,NULL,model_handle,"penalties","Matrix","Mat"); + } + + counter=0; + if(model->numpenalties){ + + /*First deal with internal grids: */ + if (model->penaltypartitioning[i]>=0){ //this penalty belongs to at least this CPU + + for(j=0;jnumlayers-1;j++){ + + /*We are pairing grids along a vertical profile.*/ + grid1=(int)*(model->penalties+model->numlayers*i+j); + grid2=(int)*(model->penalties+model->numlayers*i+j+1); + + penpair_id=counter+1; //matlab indexing + penpair_dof=1; + penpair_node_ids[0]=grid1; + penpair_node_ids[1]=grid2; + penpair_penalty_offset=model->penalty_offset; + penpair_numdofs=1; + + penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, + penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); + loads->AddObject(penpair); + + + counter++; + + penpair_id=counter+1; //matlab indexing + penpair_dof=2; + penpair_node_ids[0]=grid1; + penpair_node_ids[1]=grid2; + penpair_penalty_offset=model->penalty_offset; + penpair_numdofs=1; + + penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, + penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); + loads->AddObject(penpair); + + + counter++; + + } + } //for ( i=0; ipenalties); + + /*Create penpair loads also for rift grids, so that they don't penetrate one another, if needed: */ + /*First fetch rifts: */ + if(model->numrifts){ + + for(i=0;inumrifts;i++){ + riftpenaltypairs=model->riftspenaltypairs[i]; + for(j=0;jriftsnumpenaltypairs[i];j++){ + + el1=(int)*(riftpenaltypairs+7*j+2); + #ifdef _PARALLEL_ + if (model->epart[el1-1]!=my_rank){ + /*This load does not belong to this cluster node, as it references an element + *that does not belong to this node's partition. Drop this 'j':*/ + continue; + } + #endif + + /*Ok, retrieve all the data needed to add a penalty between the two grids: */ + el2=(int)*(riftpenaltypairs+7*j+3); + grid1=(int)*(riftpenaltypairs+7*j+0); + grid2=(int)*(riftpenaltypairs+7*j+1); + normal[0]=*(riftpenaltypairs+7*j+4); + normal[1]=*(riftpenaltypairs+7*j+5); + length=*(riftpenaltypairs+7*j+6); + friction=model->riftsfriction[i]; + fill=model->riftsfill[i]; + + penpair_id=counter+1; //matlab indexing + penpair_node_ids[0]=grid1; + penpair_node_ids[1]=grid2; + penpair_element_ids[0]=el1; + penpair_element_ids[1]=el2; + penpair_penalty_offset=model->penalty_offset; + penpair_penalty_lock=model->penalty_lock; + penpair_numdofs=2; + penpair_normal[0]=normal[0]; + penpair_normal[1]=normal[1]; + penpair_length=length; + penpair_friction=friction; + penpair_fill=fill; + + penpair= new Penpair(penpair_id,penpair_penalty_offset,penpair_penalty_lock,penpair_numdofs,penpair_node_ids,penpair_dof, + penpair_element_ids,penpair_friction,penpair_fill,penpair_normal,penpair_length); + + loads->AddObject(penpair); + + counter++; + } + } + } + + /*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(); + + /*Free ressources:*/ + xfree((void**)&riftsnumpenaltypairs); + for(i=0;inumrifts;i++){ + double* temp=riftspenaltypairs[i]; + xfree((void**)&temp); + } + xfree((void**)&riftspenaltypairs); + xfree((void**)&riftsfill); + xfree((void**)&riftsfriction); + + + /*Assign output pointer: */ + *ploads=loads; + +} + + Index: sm/trunk/src/c/ModelProcessorx/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DistributeNumDofs.cpp (revision 116) +++ (revision ) @@ -1,33 +1,0 @@ -/*!\file: DistributeNumDofs.cpp - * \brief: figure out the maximum number of dofs per grid. - */ - -#include "../include/macros.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../shared/shared.h" - -#undef __FUNCT__ -#define __FUNCT__ "DistributeNumDofs" - -int DistributeNumDofs(int *pnumdofs,int analysis_type){ - - int numdofs=2; - int i; - - /*ok, according to analysis type: */ - if (analysis_type==DiagnosticHorizAnalysisEnum()){ - numdofs=2; - } - else if (analysis_type==ControlAnalysisEnum()){ - numdofs=2; - } - else if (analysis_type==ThermalSteadyAnalysisEnum()){ - numdofs=1; - } - else throw ErrorException(__FUNCT__," analysis type not implemented yet!"); - - /*Assign output pointers:*/ - *pnumdofs=numdofs;; - - return 1; -} Index: sm/trunk/src/c/ModelProcessorx/FetchRifts.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/FetchRifts.cpp (revision 116) +++ (revision ) @@ -1,78 +1,0 @@ -/*! \file FetchRifts.cpp - * \brief: special i/o here to recover the rifts data. Because serially, the rifts are held into a matlab array strucutre, - * hard to parse. In parallel, it's easier, as we have marshalled the data into a binary buffer. - */ - -#undef __FUNCT__ -#define __FUNCT__ "FetchRifts" - -#include "../io/io.h" -#include "../shared/shared.h" - -int FetchRifts(int** priftsnumpenaltypairs,double*** priftspenaltypairs,int** priftsfill,double** priftsfriction,ConstDataHandle model_handle,int numrifts){ - - int i; - int noerr=1; - - /*output: */ - int* riftsnumpenaltypairs=NULL; - double** riftspenaltypairs=NULL; - int* riftsfill=NULL; - double* riftsfriction=NULL; - - /*intermediary: */ - double fill; - double* riftpenaltypairs=NULL; - #ifdef _SERIAL_ - mxArray* pmxa_rifts=NULL; - #endif - - if(numrifts){ - - /*Allocate arrays: */ - riftsnumpenaltypairs=(int*)xmalloc(numrifts*sizeof(int)); - riftspenaltypairs=(double**)xmalloc(numrifts*sizeof(double*)); - riftsfill=(int*)xmalloc(numrifts*sizeof(int)); - riftsfriction=(double*)xmalloc(numrifts*sizeof(double)); - - #ifdef _SERIAL_ - /*From model handle, recover rifts handle: */ - pmxa_rifts=mxGetField(model_handle,0,"rifts"); - - for(i=0;ianalysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ - - CreateElementsNodesAndMaterialsDiagnosticHoriz(pelements,pnodes,pmaterials, model,model_handle); - - } - /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){ - - noerr=CreateElementsNodesAndMaterialsDataSetsDiagnosticBaseVert(pelements,pnodes,pmaterials, model,model->analysis_type); - - } - else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ - - noerr=CreateElementsNodesAndMaterialsDataSetsDiagnosticVert(pelements,pnodes,pmaterials, model,model->analysis_type); - } - else if (strcmp(model->analysis_type,"melting")==0){ - - noerr=CreateElementsNodesAndMaterialsDataSetsMelting(pelements,pnodes,pmaterials, model,model->analysis_type); - } - else if (strcmp(model->analysis_type,"prognostic")==0){ - - noerr=CreateElementsNodesAndMaterialsDataSetsPrognostic(pelements,pnodes,pmaterials, model,model->analysis_type); - } - else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ - - noerr=CreateElementsNodesAndMaterialsDataSetsThermal(pelements,pnodes,pmaterials, model,model->analysis_type); - }*/ - else{ - _printf_("%s%s\n",__FUNCT__,"error message: analysis_type not supported yet!"); - } - - return noerr; - -} - - -/*! ModelCreateConstraints - */ - -#undef __FUNCT__ -#define __FUNCT__ "ModelCreateConstraints" - -int ModelCreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle){ - - int noerr=1; - - /*This is just a high level driver: */ - if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ - noerr=CreateConstraintsDiagnosticHoriz(pconstraints,model,model_handle); - } - /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){ - noerr=CreateConstraintsDiagnosticBaseVert(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ - noerr=CreateConstraintsDiagnosticVert(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"melting")==0){ - noerr=CreateConstraintsMelting(pconstraints,model,model_handle); - } - else if (strcmp(model->analysis_type,"prognostic")==0){ - noerr=CreateConstraintsPrognostic(pconstraints,model,model_handle); - } - else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ - noerr=CreateConstraintsThermal(pconstraints,model,model_handle); - } - else{ - _printf_("%s%s\n",__FUNCT__,"error message: model->analysis_type not supported yet!"); - }*/ - - return noerr; - -} - - - -/*! ModelCreateLoads - */ - - -#undef __FUNCT__ -#define __FUNCT__ "ModelCreateLoads" - -int ModelCreateLoads(DataSet** ploads, Model* model,ConstDataHandle model_handle){ - - int noerr=1; - - /*This is just a high level driver: */ - if ((strcmp(model->analysis_type,"diagnostic_horiz")==0)|| (strcmp(model->analysis_type,"control")==0)){ - - noerr=CreateLoadsDiagnosticHoriz(ploads,model,model_handle); - } - /*else if (strcmp(model->analysis_type,"diagnostic_basevert")==0){ - - noerr=CreateLoadsDiagnosticBaseVert(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"diagnostic_vert")==0){ - - noerr=CreateLoadsDiagnosticVert(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"melting")==0){ - - noerr=CreateLoadsMelting(ploads,model,model_handle); - } - else if (strcmp(model->analysis_type,"prognostic")==0){ - - noerr=CreateLoadsPrognostic(ploads,model,model_handle); - } - else if ((strcmp(model->analysis_type,"thermalsteady")==0) || (strcmp(model->analysis_type,"thermaltransient")==0)){ - - noerr=CreateLoadsThermal(ploads,model,model_handle); - - } - else{ - _printf_("%s%s\n",__FUNCT__,"error message: model->analysis_type not supported yet!"); - }*/ - - return noerr; - -} Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 117) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 117) @@ -0,0 +1,181 @@ +/* \file Mode.h + * \brief Header file defining the Model structure and processing of input data. + * \sa Model.cpp + */ + +#ifndef _MODEL_H +#define _MODEL_H + +#include "../io/io.h" +#include "../DataSet/DataSet.h" +#include "../toolkits/toolkits.h" + +struct Model { + + char* repository; + char* meshtype; + char* analysis_type; + char* solverstring; + + /*2d mesh: */ + int numberofelements; + int numberofnodes; + double* x; + double* y; + double* z; + double* elements; + double* elements_type; + + /*3d mesh: */ + int numberofnodes2d; + int numberofelements2d; + double* elements2d; + double* deadgrids; + int numlayers; + double* uppernodes; + + + /*results: */ + double* vx; + double* vy; + double* vz; + double* pressure; + + /*observations: */ + double* vx_obs; + double* vy_obs; + + /*geometry: */ + double* elementonbed; + double* elementonsurface; + double* gridonbed; + double* gridonsurface; + double* thickness; + double* surface; + double* bed; + double* elementoniceshelf; + + /*friction: */ + int drag_type; + double* drag; + double* p; + double* q; + + /*boundary conditions: */ + //diagnostic + int numberofsegs_diag; + double* segmentonneumann_diag; + double* neumannvalues_diag; + double* gridondirichlet_diag; + double* dirichletvalues_diag; + //prognostic + double* segmentonneumann_prog; + double* neumannvalues_prog; + double* gridondirichlet_prog; + double* dirichletvalues_prog; + //prognostic2 + double* segmentonneumann_prog2; + double* neumannvalues_prog2; + //thermal + double* gridondirichlet_thermal; + double* dirichletvalues_thermal; + double* geothermalflux; + + /*materials: */ + double rho_water,rho_ice; + double g; + double* B; + double* n; + + /*control methods: */ + char* control_type; + + /*solution parameters: */ + double* fit; + double meanvel,epsvel; + int artificial_diffusivity; + int nsteps; + double tolx; + double* maxiter; + double mincontrolconstraint; + double maxcontrolconstraint; + int debug; + int plot; + double eps_rel; + double eps_abs; + double dt,ndt; + double penalty_offset; + double penalty_melting; + int penalty_lock; + double sparsity; + int connectivity; + int lowmem; + double* optscal; + double yts; + double viscosity_overshoot; + int waitonlock; + + /*thermal parameters: */ + double beta; + double meltingpoint; + double latentheat; + double heatcapacity,thermalconductivity; + int min_thermal_constraints; + double mixed_layer_capacity; + double thermal_exchange_velocity; + + /*rifts: */ + int numrifts; + int* riftsnumpenaltypairs; + double** riftspenaltypairs; + int* riftsfill; + double* riftsfriction; + + /*penalties: */ + int numpenalties; + double* penalties; + + /*basal: */ + double* melting; + double* accumulation; + + + /*exterior data: */ + int* epart; /*!element partitioning.*/ + int* my_grids; /*! grids that belong to this cpu*/ + double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/ + int* penaltypartitioning; + +}; + + + + /*constructor and destructor: */ + Model* NewModel(void); + void DeleteModel( Model** pthis); + + /*Echo: */ + void ModelEcho(Model* model,int which_part,int rank); + + /*Initialization with matlab workspace data, or marshall binary data: */ + int ModelInit(Model** pmodel,ConstDataHandle model_handle); + + /*Creation of fem datasets: general drivers*/ + void CreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); + void CreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle); + void CreateLoads(DataSet** ploads, Model* model, ConstDataHandle model_handle); + void CreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle); + + + /*Create of fem datasets: specialised drivers: */ + + /*diagnostic:*/ + void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); + void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); + void CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle); + + /*control:*/ + void CreateParametersControl(DataSet* parameters,Model* model,ConstDataHandle model_handle,int* pcount); + + +#endif /* _MODEL_H */ Index: sm/trunk/src/c/ModelProcessorx/ModelCreateParameters.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/ModelCreateParameters.cpp (revision 116) +++ (revision ) @@ -1,280 +1,0 @@ -/* - * ModelCreateParameters.c: - */ - -#undef __FUNCT__ -#define __FUNCT__ "ModelCreateParameters" - -#include "../DataSet/DataSet.h" -#include "../toolkits/toolkits.h" -#include "../EnumDefinitions/EnumDefinitions.h" -#include "../objects/objects.h" -#include "../shared/shared.h" -#include "./ModelProcessorx.h" - -int ModelCreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle){ - - int noerr=1; - int i; - - DataSet* parameters = NULL; - Param* param = NULL; - int count=1; - int analysis_type; - int numberofdofspernode; - - double* fit=NULL; - double* optscal=NULL; - double* maxiter=NULL; - double* control_parameter=NULL; - - - - /*Initialize dataset: */ - parameters = new DataSet(ParametersEnum()); - - //Get analysis_type: - analysis_type=AnalysisTypeAsEnum(model->analysis_type); - - count++; - param= new Param(count,"analysis_type",INTEGER); - param->SetInteger(analysis_type); - parameters->AddObject(param); - - - /*debug: */ - param= new Param(count,"debug",INTEGER); - param->SetInteger(model->debug); - parameters->AddObject(param); - - /*eps_rel: */ - count++; - param= new Param(count,"eps_rel",DOUBLE); - param->SetDouble(model->eps_rel); - parameters->AddObject(param); - - /*eps_abs: */ - count++; - param= new Param(count,"eps_abs",DOUBLE); - param->SetDouble(model->eps_abs); - parameters->AddObject(param); - - /*yts: */ - count++; - param= new Param(count,"yts",DOUBLE); - param->SetDouble(model->yts); - parameters->AddObject(param); - - /*dt: */ - count++; - param= new Param(count,"dt",DOUBLE); - param->SetDouble(model->dt); - parameters->AddObject(param); - - /*ndt: */ - count++; - param= new Param(count,"ndt",DOUBLE); - param->SetDouble(model->ndt); - parameters->AddObject(param); - - /*penalty_offset: */ - count++; - param= new Param(count,"penalty_offset",DOUBLE); - param->SetDouble(model->penalty_offset); - parameters->AddObject(param); - - /*sparsity: */ - count++; - param= new Param(count,"sparsity",DOUBLE); - param->SetDouble(model->sparsity); - parameters->AddObject(param); - - /*lowmem: */ - count++; - param= new Param(count,"lowmem",INTEGER); - param->SetInteger(model->lowmem); - parameters->AddObject(param); - - /*connectivity: */ - count++; - param= new Param(count,"connectivity",INTEGER); - param->SetInteger(model->connectivity); - parameters->AddObject(param); - - /*beta: */ - count++; - param= new Param(count,"beta",DOUBLE); - param->SetDouble(model->beta); - parameters->AddObject(param); - - /*meltingpoint: */ - count++; - param= new Param(count,"meltingpoint",DOUBLE); - param->SetDouble(model->meltingpoint); - parameters->AddObject(param); - - /*latentheat: */ - count++; - param= new Param(count,"latentheat",DOUBLE); - param->SetDouble(model->latentheat); - parameters->AddObject(param); - - /*heatcapacity: */ - count++; - param= new Param(count,"heatcapacity",DOUBLE); - param->SetDouble(model->heatcapacity); - parameters->AddObject(param); - - /*penalty_melting: */ - count++; - param= new Param(count,"penalty_melting",DOUBLE); - param->SetDouble(model->penalty_melting); - parameters->AddObject(param); - - /*min_thermal_constraints: */ - count++; - param= new Param(count,"min_thermal_constraints",INTEGER); - param->SetInteger(model->min_thermal_constraints); - parameters->AddObject(param); - - /*waitonlock: */ - count++; - param= new Param(count,"waitonlock",INTEGER); - param->SetInteger(model->waitonlock); - parameters->AddObject(param); - - /*solverstring: */ - count++; - param= new Param(count,"solverstring",STRING); - param->SetString(model->solverstring); - parameters->AddObject(param); - - /*plot: */ - count++; - param= new Param(count,"plot",INTEGER); - param->SetInteger(model->plot); - parameters->AddObject(param); - - /*numberofgrids: */ - count++; - param= new Param(count,"numberofnodes",INTEGER); - param->SetInteger(model->numberofnodes); - parameters->AddObject(param); - - /*Deal with numberofdofspernode: */ - DistributeNumDofs(&numberofdofspernode,analysis_type); - - count++; - param= new Param(count,"numberofdofspernode",INTEGER); - param->SetInteger(numberofdofspernode); - parameters->AddObject(param); - - /*All our datasets are already ordered by ids. Set presort flag so that later on, when sorting is requested on these - * datasets, it will not be redone: */ - parameters->Presort(); - - /*Deal with control parameters: */ - if (analysis_type==ControlAnalysisEnum()){ - - /*control_type: */ - count++; - param= new Param(count,"control_type",STRING); - param->SetString(model->control_type); - parameters->AddObject(param); - - - /*nsteps: */ - count++; - param= new Param(count,"nsteps",INTEGER); - param->SetInteger(model->nsteps); - parameters->AddObject(param); - - /*tolx: */ - count++; - param= new Param(count,"tolx",DOUBLE); - param->SetDouble(model->tolx); - parameters->AddObject(param); - - /*mincontrolconstraint: */ - count++; - param= new Param(count,"mincontrolconstraint",DOUBLE); - param->SetDouble(model->mincontrolconstraint); - parameters->AddObject(param); - - /*maxcontrolconstraint: */ - count++; - param= new Param(count,"maxcontrolconstraint",DOUBLE); - param->SetDouble(model->maxcontrolconstraint); - parameters->AddObject(param); - - /*epsvel: */ - count++; - param= new Param(count,"epsvel",DOUBLE); - param->SetDouble(model->epsvel); - parameters->AddObject(param); - - /*meanvel: */ - count++; - param= new Param(count,"meanvel",DOUBLE); - param->SetDouble(model->meanvel); - parameters->AddObject(param); - - /*Now, recover fit, optscal and maxiter as vectors: */ - ModelFetchData((void**)&model->fit,NULL,NULL,model_handle,"fit","Matrix","Mat"); - ModelFetchData((void**)&model->optscal,NULL,NULL,model_handle,"optscal","Matrix","Mat"); - ModelFetchData((void**)&model->maxiter,NULL,NULL,model_handle,"maxiter","Matrix","Mat"); - - count++; - param= new Param(count,"fit",DOUBLEVEC); - param->SetDoubleVec(model->fit,model->nsteps); - parameters->AddObject(param); - - count++; - param= new Param(count,"optscal",DOUBLEVEC); - param->SetDoubleVec(model->optscal,model->nsteps); - parameters->AddObject(param); - - count++; - param= new Param(count,"maxiter",DOUBLEVEC); - param->SetDoubleVec(model->maxiter,model->nsteps); - parameters->AddObject(param); - - xfree((void**)&model->fit); - xfree((void**)&model->optscal); - xfree((void**)&model->maxiter); - - /*Get vx_obs, vy_obs, and the parameter value: */ - ModelFetchData((void**)&model->vx_obs,NULL,NULL,model_handle,"vx_obs","Matrix","Mat"); - ModelFetchData((void**)&model->vy_obs,NULL,NULL,model_handle,"vy_obs","Matrix","Mat"); - ModelFetchData((void**)&control_parameter,NULL,NULL,model_handle,model->control_type,"Matrix","Mat"); - - for(i=0;inumberofnodes;i++)model->vx_obs[i]=model->vx_obs[i]/model->yts; - for(i=0;inumberofnodes;i++)model->vy_obs[i]=model->vy_obs[i]/model->yts; - - count++; - param= new Param(count,"vx_obs",DOUBLEVEC); - param->SetDoubleVec(model->vx_obs,model->numberofnodes); - parameters->AddObject(param); - - count++; - param= new Param(count,"vy_obs",DOUBLEVEC); - param->SetDoubleVec(model->vy_obs,model->numberofnodes); - parameters->AddObject(param); - - count++; - param= new Param(count,"control_parameter",DOUBLEVEC); - param->SetDoubleVec(control_parameter,model->numberofnodes); - parameters->AddObject(param); - - xfree((void**)&model->vx_obs); - xfree((void**)&model->vy_obs); - xfree((void**)&control_parameter); - - - } - - /*Assign output pointer: */ - *pparameters=parameters; -} - - Index: sm/trunk/src/c/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/ModelProcessorx.h (revision 116) +++ (revision ) @@ -1,215 +1,0 @@ -/* \file ModelProcessorx.h - * \brief Header file defining the Model structure and processing of input data. - * \sa Model.cpp - */ - -#ifndef _MODELPROCESSORX_H -#define _MODELPROCESSORX_H - -#include "../io/io.h" -#include "../DataSet/DataSet.h" -#include "../toolkits/toolkits.h" - -typedef struct { - - char* repository; - char* meshtype; - char* analysis_type; - char* solverstring; - - /*2d mesh: */ - int numberofelements; - int numberofnodes; - double* x; - double* y; - double* z; - double* elements; - double* elements_type; - - /*3d mesh: */ - int numberofnodes2d; - int numberofelements2d; - double* elements2d; - double* deadgrids; - int numlayers; - double* uppernodes; - - - /*results: */ - double* vx; - double* vy; - double* vz; - double* pressure; - - /*observations: */ - double* vx_obs; - double* vy_obs; - - /*geometry: */ - double* elementonbed; - double* elementonsurface; - double* gridonbed; - double* gridonsurface; - double* thickness; - double* surface; - double* bed; - double* elementoniceshelf; - - /*friction: */ - int drag_type; - double* drag; - double* p; - double* q; - - /*boundary conditions: */ - //diagnostic - int numberofsegs_diag; - double* segmentonneumann_diag; - double* neumannvalues_diag; - double* gridondirichlet_diag; - double* dirichletvalues_diag; - //prognostic - double* segmentonneumann_prog; - double* neumannvalues_prog; - double* gridondirichlet_prog; - double* dirichletvalues_prog; - //prognostic2 - double* segmentonneumann_prog2; - double* neumannvalues_prog2; - //thermal - double* gridondirichlet_thermal; - double* dirichletvalues_thermal; - double* geothermalflux; - - /*materials: */ - double rho_water,rho_ice; - double g; - double* B; - double* n; - - /*control methods: */ - char* control_type; - - /*solution parameters: */ - double* fit; - double meanvel,epsvel; - int artificial_diffusivity; - int nsteps; - double tolx; - double* maxiter; - double mincontrolconstraint; - double maxcontrolconstraint; - int debug; - int plot; - double eps_rel; - double eps_abs; - double dt,ndt; - double penalty_offset; - double penalty_melting; - int penalty_lock; - double sparsity; - int connectivity; - int lowmem; - double* optscal; - double yts; - double viscosity_overshoot; - int waitonlock; - - /*thermal parameters: */ - double beta; - double meltingpoint; - double latentheat; - double heatcapacity,thermalconductivity; - int min_thermal_constraints; - double mixed_layer_capacity; - double thermal_exchange_velocity; - - /*rifts: */ - int numrifts; - int* riftsnumpenaltypairs; - double** riftspenaltypairs; - int* riftsfill; - double* riftsfriction; - - /*penalties: */ - int numpenalties; - double* penalties; - - /*basal: */ - double* melting; - double* accumulation; - - - /*exterior data: */ - int* epart; /*!element partitioning.*/ - int* my_grids; /*! grids that belong to this cpu*/ - double* my_bordergrids; /*! grids that belong to this cpu, and some other cpu also*/ - int* penaltypartitioning; - -} Model; - - - - Model* NewModel(void); - - void DeleteModel( Model** pthis); - - void ModelEcho(Model* model,int which_part,int rank); - - int ModelInit(Model** pmodel,ConstDataHandle model_handle); - - int ModelCreateElementsNodesAndMaterials(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); - int CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, Model* model,ConstDataHandle model_handle); - int DistributeNumDofs(int *pnumdofs,int analysis_type); - - int ModelCreateConstraints(DataSet** pconstraints, Model* model,ConstDataHandle model_handle); - int CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,Model* model,ConstDataHandle model_handle); - - int ModelCreateLoads(DataSet** ploads, Model* model, ConstDataHandle model_handle); - int CreateLoadsDiagnosticHoriz(DataSet** ploads, Model* model, ConstDataHandle model_handle); - - int ModelCreateParameters(DataSet** pparameters,Model* model,ConstDataHandle model_handle); - -#ifdef _PARALLEL_ - //int ProcessBatchParameters(BatchParams* params,Model* model,DataSet* ncase); - //int ProcessWorkspaceParameters(WorkspaceParams* params,Model* model,Vec* partition); -#endif - -int FetchRifts(int** priftsnumpenaltypairs,double*** priftspenaltypairs,int** priftsfill,double** priftsfriction,ConstDataHandle model_handle,int numrifts); - -/*Elements, grids and materials: */ -int CreateElementsGridsAndMaterialsDataSets(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); -int CreateElementsGridsAndMaterialsDataSetsDiagnosticBaseVert(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); -int CreateElementsGridsAndMaterialsDataSetsDiagnosticHoriz(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); -int CreateElementsGridsAndMaterialsDataSetsDiagnosticVert(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); -int CreateElementsGridsAndMaterialsDataSetsMelting(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); -int CreateElementsGridsAndMaterialsDataSetsPrognostic(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); -int CreateElementsGridsAndMaterialsDataSetsThermal(int** pepart, DataSet** pelements, DataSet** pgrids, DataSet** pmaterials, Vec* ppartition, Vec* ptpartition, - int** pmy_grids, Vec* pmy_bordergrids, Model* model,char* analysis_type); - -/*Constraints: */ -int CreateConstraintsDataSet(DataSet** pconstraints,Model* model,char* analysis_type); -int CreateConstraintsDataSetDiagnosticBaseVert(DataSet** pconstraints,Model* model,char* analysis_type); -int CreateConstraintsDataSetDiagnosticHoriz(DataSet** pconstraints,Model* model,char* analysis_type); -int CreateConstraintsDataSetDiagnosticVert(DataSet** pconstraints,Model* model,char* analysis_type); -int CreateConstraintsDataSetMelting(DataSet** pconstraints,Model* model,char* analysis_type); -int CreateConstraintsDataSetPrognostic(DataSet** pconstraints,Model* model,char* analysis_type); -int CreateConstraintsDataSetThermal(DataSet** pconstraints,Model* model,char* analysis_type); - -/*Loads: */ -int CreateLoadsDataSet(DataSet** ploads,int* my_grids,double* my_bordergrids, int* epart,Model* model,char* analysis_type); -int CreateLoadsDataSetDiagnosticBaseVert(DataSet** ploads,int* my_grids,double* my_bordergrids,int* epart,Model* model,char* analysis_type); -int CreateLoadsDataSetDiagnosticHoriz(DataSet** ploads,int* my_grids,double* my_bordergrids,int* epart,Model* model,char* analysis_type); -int CreateLoadsDataSetDiagnosticVert(DataSet** ploads,int* my_grids,double* my_bordergrids,int* epart,Model* model,char* analysis_type); -int CreateLoadsDataSetMelting(DataSet** ploads,int* my_grids,double* my_bordergrids,int* epart,Model* model,char* analysis_type); -int CreateLoadsDataSetPrognostic(DataSet** ploads,int* my_grids,double* my_bordergrids,int* epart,Model* model,char* analysis_type); -int CreateLoadsDataSetThermal(DataSet** ploads,int* my_grids,double* my_bordergrids,int* epart,Model* model,char* analysis_type); - - -#endif /* _MODELPROCESSORX_H */ Index: /issm/trunk/src/c/io/FetchRifts.cpp =================================================================== --- /issm/trunk/src/c/io/FetchRifts.cpp (revision 117) +++ /issm/trunk/src/c/io/FetchRifts.cpp (revision 117) @@ -0,0 +1,78 @@ +/*! \file FetchRifts.cpp + * \brief: special i/o here to recover the rifts data. Because serially, the rifts are held into a matlab array strucutre, + * hard to parse. In parallel, it's easier, as we have marshalled the data into a binary buffer. + */ + +#undef __FUNCT__ +#define __FUNCT__ "FetchRifts" + +#include "./io.h" +#include "../shared/shared.h" + +int FetchRifts(int** priftsnumpenaltypairs,double*** priftspenaltypairs,int** priftsfill,double** priftsfriction,ConstDataHandle model_handle,int numrifts){ + + int i; + int noerr=1; + + /*output: */ + int* riftsnumpenaltypairs=NULL; + double** riftspenaltypairs=NULL; + int* riftsfill=NULL; + double* riftsfriction=NULL; + + /*intermediary: */ + double fill; + double* riftpenaltypairs=NULL; + #ifdef _SERIAL_ + mxArray* pmxa_rifts=NULL; + #endif + + if(numrifts){ + + /*Allocate arrays: */ + riftsnumpenaltypairs=(int*)xmalloc(numrifts*sizeof(int)); + riftspenaltypairs=(double**)xmalloc(numrifts*sizeof(double*)); + riftsfill=(int*)xmalloc(numrifts*sizeof(int)); + riftsfriction=(double*)xmalloc(numrifts*sizeof(double)); + + #ifdef _SERIAL_ + /*From model handle, recover rifts handle: */ + pmxa_rifts=mxGetField(model_handle,0,"rifts"); + + for(i=0;i