Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 2721) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.cpp (revision 2722) @@ -42,4 +42,5 @@ int PrognosticAnalysisEnum(void){ return 250; } int BalancedthicknessAnalysisEnum(void){ return 251; } +int BalancedvelocitiesAnalysisEnum(void){ return 252; } //melting int MeltingAnalysisEnum(void){ return 260; } Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 2721) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 2722) @@ -43,4 +43,5 @@ int PrognosticAnalysisEnum(void); int BalancedthicknessAnalysisEnum(void); +int BalancedvelocitiesAnalysisEnum(void); //melting int MeltingAnalysisEnum(void); Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 2721) +++ /issm/trunk/src/c/Makefile.am (revision 2722) @@ -222,4 +222,8 @@ ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\ ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp\ ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\ ./Dofx/Dofx.h\ @@ -525,4 +529,8 @@ ./ModelProcessorx/Balancedthickness/CreateLoadsBalancedthickness.cpp\ ./ModelProcessorx/Balancedthickness/CreateParametersBalancedthickness.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp\ + ./ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp\ ./ModelProcessorx/Qmu/CreateParametersQmu.cpp\ ./Dofx/Dofx.h\ @@ -625,4 +633,5 @@ ./parallel/prognostic_core.cpp\ ./parallel/balancedthickness_core.cpp\ + ./parallel/balancedvelocities_core.cpp\ ./parallel/transient_core.cpp\ ./parallel/transient_core_2d.cpp\ @@ -636,5 +645,5 @@ bin_PROGRAMS = else -bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe transient.exe steadystate.exe +bin_PROGRAMS = diagnostic.exe thermal.exe prognostic.exe balancedthickness.exe balancedvelocities.exe transient.exe steadystate.exe endif @@ -656,4 +665,7 @@ balancedthickness_exe_CXXFLAGS= -fPIC -D_PARALLEL_ +balancedvelocities_exe_SOURCES = parallel/balancedvelocities.cpp +balancedvelocities_exe_CXXFLAGS= -fPIC -D_PARALLEL_ + transient_exe_SOURCES = parallel/transient.cpp transient_exe_CXXFLAGS= -fPIC -D_PARALLEL_ Index: /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp (revision 2722) +++ /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateConstraintsBalancedvelocities.cpp (revision 2722) @@ -0,0 +1,77 @@ +/* + * CreateConstraintsBalancedvelocities.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateConstraintsBalancedvelocities" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../IoModel.h" + + +void CreateConstraintsBalancedvelocities(DataSet** pconstraints, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + + int i; + int count; + + DataSet* constraints = NULL; + Spc* spc = NULL; + + /*spc intermediary data: */ + int spc_sid; + int spc_node; + int spc_dof; + double spc_value; + + double* spcvelocity=NULL; + + /*Create constraints: */ + constraints = new DataSet(ConstraintsEnum()); + + /*Fetch data: */ + IoModelFetchData(&spcvelocity,NULL,NULL,iomodel_handle,"spcvelocity"); + + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((iomodel->my_grids[i]==1)){ + #endif + + if ((int)spcvelocity[6*i+0] && (int)spcvelocity[6*i+1]){ //spc if vx and vy are constrained + + /*This grid needs to be spc'd: */ + spc_sid=count; + spc_node=i+1; + spc_dof=1; //we enforce first translation degree of freedom, for temperature + spc_value=pow( pow(*(spcvelocity+6*i+4),2.0) + pow(*(spcvelocity+6*i+5),2.0) ,0.5); + + 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**)&spcvelocity); + + cleanup_and_return: + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (revision 2722) +++ /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateElementsNodesAndMaterialsBalancedvelocities.cpp (revision 2722) @@ -0,0 +1,541 @@ +/* + * CreateElementsNodesAndMaterialsBalancedvelocities.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateElementsNodesAndMaterialsBalancedvelocities" + +#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 "../IoModel.h" + + +void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_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; + + /*output: */ + int* epart=NULL; //element partitioning. + int* npart=NULL; //node partitioning. + int* my_grids=NULL; + double* my_bordergrids=NULL; + + + /*intermediary: */ + int elements_width; //size of elements + double B_avg; + + /*tria constructor input: */ + int tria_id; + int tria_mid; + int tria_mparid; + int tria_numparid; + int tria_g[3]; + double tria_h[3]; + double tria_s[3]; + double tria_b[3]; + double tria_k[3]; + double tria_melting[3]; + double tria_accumulation[3]; + double tria_geothermalflux[3]; + int tria_friction_type; + double tria_p; + double tria_q; + int tria_shelf; + bool tria_onwater; + + /*matice constructor input: */ + int matice_mid; + double matice_B; + double matice_n; + + /*penta constructor input: */ + + int penta_id; + int penta_mid; + int penta_mparid; + int penta_numparid; + int penta_g[6]; + double penta_h[6]; + double penta_s[6]; + double penta_b[6]; + double penta_k[6]; + int penta_friction_type; + double penta_p; + double penta_q; + int penta_shelf; + int penta_onbed; + int penta_onsurface; + int penta_collapse; + double penta_melting[6]; + double penta_accumulation[6]; + double penta_geothermalflux[6]; + int penta_thermal_steadystate; + bool penta_onwater; + + /*matpar constructor input: */ + int matpar_mid; + double matpar_rho_ice; + double matpar_rho_water; + double matpar_heatcapacity; + double matpar_thermalconductivity; + double matpar_latentheat; + double matpar_beta; + double matpar_meltingpoint; + double matpar_mixed_layer_capacity; + double matpar_thermal_exchange_velocity; + double matpar_g; + + /* node constructor input: */ + int node_id; + int node_partitionborder=0; + double node_x[3]; + double node_sigma; + int node_onbed; + int node_onsurface; + int node_onshelf; + int node_onsheet; + int node_upper_node_id; + int node_numdofs; + + + #ifdef _PARALLEL_ + /*Metis partitioning: */ + int range; + Vec gridborder=NULL; + int my_numgrids; + int* all_numgrids=NULL; + int gridcount; + int count; + #endif + int first_grid_index; + + /*First create the elements, nodes and material properties: */ + elements = new DataSet(ElementsEnum()); + nodes = new DataSet(NodesEnum()); + materials = new DataSet(MaterialsEnum()); + + /*Width of elements: */ + if(strcmp(iomodel->meshtype,"2d")==0){ + elements_width=3; //tria elements + } + else{ + elements_width=6; //penta elements + } + + #ifdef _PARALLEL_ + /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/ + if(strcmp(iomodel->meshtype,"2d")==0){ + /*load elements: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + } + else{ + /*load elements2d: */ + IoModelFetchData(&iomodel->elements2d,NULL,NULL,iomodel_handle,"elements2d"); + } + + + MeshPartitionx(&epart, &npart,iomodel->numberofelements,iomodel->numberofnodes,iomodel->elements, iomodel->numberofelements2d,iomodel->numberofnodes2d,iomodel->elements2d,iomodel->numlayers,elements_width, iomodel->meshtype,num_procs); + + /*Free elements and elements2d: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->elements2d); + + /*Used later on: */ + my_grids=(int*)xcalloc(iomodel->numberofnodes,sizeof(int)); + #endif + + + + /*elements created vary if we are dealing with a 2d mesh, or a 3d mesh: */ + + /*2d mesh: */ + if (strcmp(iomodel->meshtype,"2d")==0){ + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + + for (i=0;inumberofelements;i++){ + + #ifdef _PARALLEL_ + /*!All elements have been partitioned above, only create elements for this CPU: */ + if(my_rank==epart[i]){ + #endif + + + /*ids: */ + tria_id=i+1; //matlab indexing. + tria_mid=-1; //no need for materials + tria_mparid=-1; //no need for materials + tria_numparid=1; + + /*vertices offsets: */ + tria_g[0]=(int)*(iomodel->elements+elements_width*i+0); + tria_g[1]=(int)*(iomodel->elements+elements_width*i+1); + tria_g[2]=(int)*(iomodel->elements+elements_width*i+2); + + /*thickness,surface and bed:*/ + tria_h[0]= *(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+0)-1)); //remember, elements is an index of vertices offsets, in matlab indexing. + tria_h[1]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+1)-1)); + tria_h[2]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+2)-1)) ; + + tria_s[0]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+0)-1)); + tria_s[1]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+1)-1)); + tria_s[2]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+2)-1)); + + tria_b[0]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+0)-1)); + tria_b[1]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+1)-1)); + tria_b[2]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+2)-1)); + + /*element on iceshelf?:*/ + tria_shelf=(int)*(iomodel->elementoniceshelf+i); + tria_onwater=(bool)*(iomodel->elementonwater+i); + + /*Create tria element using its constructor:*/ + tria=new Tria(tria_id, tria_mid, tria_mparid, tria_numparid,tria_g, tria_h, tria_s, tria_b, tria_k, tria_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_onwater); + + /*Add tria element to elements dataset: */ + elements->AddObject(tria); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonwater); + + } + else{ // if (strcmp(meshtype,"2d")==0) + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + + for (i=0;inumberofelements;i++){ + #ifdef _PARALLEL_ + /*We are using our element partition to decide which elements will be created on this node: */ + if(my_rank==epart[i]){ + #endif + + + /*name and id: */ + penta_id=i+1; //matlab indexing. + penta_mid=-1; + penta_mparid=-1; //no need for materials + penta_numparid=1; + + /*vertices,thickness,surface,bed and drag: */ + for(j=0;j<6;j++){ + penta_g[j]=(int)*(iomodel->elements+elements_width*i+j); + penta_h[j]=*(iomodel->thickness+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + penta_s[j]=*(iomodel->surface+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + penta_b[j]=*(iomodel->bed+ ((int)*(iomodel->elements+elements_width*i+j)-1)); + } + + /*diverse: */ + penta_shelf=(int)*(iomodel->elementoniceshelf+i); + penta_onbed=(int)*(iomodel->elementonbed+i); + penta_onsurface=(int)*(iomodel->elementonsurface+i); + penta_collapse=1; + penta_onwater=(bool)*(iomodel->elementonwater+i); + + + /*Create Penta using its constructor:*/ + penta= new Penta( penta_id,penta_mid,penta_mparid,penta_numparid,penta_g,penta_h,penta_s,penta_b,penta_k,penta_friction_type, + penta_p,penta_q,penta_shelf,penta_onbed,penta_onsurface, + penta_collapse,penta_melting,penta_accumulation,penta_geothermalflux, + penta_thermal_steadystate,penta_onwater); + + /*Add penta element to elements dataset: */ + elements->AddObject(penta); + + #ifdef _PARALLEL_ + /*Now that we are here, we can also start building the list of grids belonging to this node partition: we use + *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing) + into the grid coordinates. If we start plugging 1 into my_grids for each index[n][i] (i=0:2), then my_grids + will hold which grids belong to this partition*/ + my_grids[(int)*(iomodel->elements+elements_width*i+0)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+1)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+2)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+3)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+4)-1]=1; + my_grids[(int)*(iomodel->elements+elements_width*i+5)-1]=1; + #endif + + #ifdef _PARALLEL_ + }//if(my_rank==epart[i]) + #endif + + }//for (i=0;ielements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elementonwater); + + } //if (strcmp(meshtype,"2d")==0) + + #ifdef _PARALLEL_ + /*From the element partitioning, we can determine which grids are on the inside of this cpu's + *element partition, and which are on its border with other nodes:*/ + gridborder=NewVec(iomodel->numberofnodes); + + for (i=0;inumberofnodes;i++){ + if(my_grids[i])VecSetValue(gridborder,i,1,ADD_VALUES); + } + VecAssemblyBegin(gridborder); + VecAssemblyEnd(gridborder); + + #ifdef _ISSM_DEBUG_ + VecView(gridborder,PETSC_VIEWER_STDOUT_WORLD); + #endif + + VecToMPISerial(&my_bordergrids,gridborder); + + #ifdef _ISSM_DEBUG_ + if(my_rank==0){ + for (i=0;inumberofnodes;i++){ + printf("Grid id %i Border grid %lf\n",i+1,my_bordergrids[i]); + } + } + #endif + #endif + + /*Add one constant material property to materials: */ + matpar_mid=1; //put it at the end of the materials + matpar_g=iomodel->g; + matpar_rho_ice=iomodel->rho_ice; + matpar_rho_water=iomodel->rho_water; + matpar_thermalconductivity=iomodel->thermalconductivity; + matpar_heatcapacity=iomodel->heatcapacity; + matpar_latentheat=iomodel->latentheat; + matpar_beta=iomodel->beta; + matpar_meltingpoint=iomodel->meltingpoint; + matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; + matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; + + /*Create matpar object using its constructor: */ + matpar=new Matpar(matpar_mid,matpar_rho_ice,matpar_rho_water,matpar_heatcapacity,matpar_thermalconductivity, + matpar_latentheat,matpar_beta,matpar_meltingpoint,matpar_mixed_layer_capacity, + matpar_thermal_exchange_velocity,matpar_g); + + /*Add to materials datset: */ + materials->AddObject(matpar); + + /*Partition penalties in 3d: */ + if(strcmp(iomodel->meshtype,"3d")==0){ + + /*Get penalties: */ + IoModelFetchData(&iomodel->penalties,&iomodel->numpenalties,NULL,iomodel_handle,"penalties"); + + if(iomodel->numpenalties){ + + iomodel->penaltypartitioning=(int*)xmalloc(iomodel->numpenalties*sizeof(int)); + #ifdef _SERIAL_ + for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=1; + #else + for(i=0;inumpenalties;i++)iomodel->penaltypartitioning[i]=-1; + + for(i=0;inumpenalties;i++){ + first_grid_index=(int)(*(iomodel->penalties+i*iomodel->numlayers+0)-1); + if((my_grids[first_grid_index]==1) && (my_bordergrids[first_grid_index]<=1.0) ) { //this grid belongs to this node's internal partition grids + /*All grids that are being penalised belong to this node's internal grid partition.:*/ + iomodel->penaltypartitioning[i]=1; + } + if(my_bordergrids[first_grid_index]>1.0) { //this grid belongs to a partition border + iomodel->penaltypartitioning[i]=0; + } + } + #endif + } + + /*Free penalties: */ + xfree((void**)&iomodel->penalties); + } + + /*Ok, let's summarise. Now, every CPU has the following two arrays: my_grids, and my_bordergrids. + We can therefore determine which grids are internal to this node's partition + and which ones are shared with other nodes because they are on the border of this node's partition. Knowing + that, go and create the grids*/ + + /*Create nodes from x,y,z, as well as the spc values on those grids: */ + + /*First fetch data: */ + if (strcmp(iomodel->meshtype,"3d")==0){ + IoModelFetchData(&iomodel->deadgrids,NULL,NULL,iomodel_handle,"deadgrids"); + IoModelFetchData(&iomodel->uppernodes,NULL,NULL,iomodel_handle,"uppergrids"); + } + IoModelFetchData(&iomodel->x,NULL,NULL,iomodel_handle,"x"); + IoModelFetchData(&iomodel->y,NULL,NULL,iomodel_handle,"y"); + IoModelFetchData(&iomodel->z,NULL,NULL,iomodel_handle,"z"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->gridonbed,NULL,NULL,iomodel_handle,"gridonbed"); + IoModelFetchData(&iomodel->gridonsurface,NULL,NULL,iomodel_handle,"gridonsurface"); + IoModelFetchData(&iomodel->gridonicesheet,NULL,NULL,iomodel_handle,"gridonicesheet"); + IoModelFetchData(&iomodel->gridoniceshelf,NULL,NULL,iomodel_handle,"gridoniceshelf"); + + + /*Get number of dofs per node: */ + DistributeNumDofs(&node_numdofs,iomodel->analysis_type,iomodel->sub_analysis_type); + + for (i=0;inumberofnodes;i++){ + #ifdef _PARALLEL_ + /*keep only this partition's nodes:*/ + if((my_grids[i]==1)){ + #endif + + node_id=i+1; //matlab indexing + + + + #ifdef _PARALLEL_ + if(my_bordergrids[i]>1.0) { //this grid belongs to a partition border + node_partitionborder=1; + } + else{ + node_partitionborder=0; + } + #else + node_partitionborder=0; + #endif + + node_x[0]=iomodel->x[i]; + node_x[1]=iomodel->y[i]; + node_x[2]=iomodel->z[i]; + node_sigma=(iomodel->z[i]-iomodel->bed[i])/(iomodel->thickness[i]); + + node_onbed=(int)iomodel->gridonbed[i]; + node_onsurface=(int)iomodel->gridonsurface[i]; + node_onshelf=(int)iomodel->gridoniceshelf[i]; + node_onsheet=(int)iomodel->gridonicesheet[i]; + + if (strcmp(iomodel->meshtype,"3d")==0){ + if (isnan(iomodel->uppernodes[i])){ + node_upper_node_id=node_id; //nodes on surface do not have upper nodes, only themselves. + } + else{ + node_upper_node_id=(int)iomodel->uppernodes[i]; + } + } + else{ + /*If we are running 2d, upper_node does not mean much. Just point towards itself!:*/ + node_upper_node_id=node_id; + } + + /*Create node using its constructor: */ + node=new Node(node_id,node_partitionborder,node_numdofs,node_x,node_sigma,node_onbed,node_onsurface,node_upper_node_id,node_onshelf,node_onsheet); + + /*set single point constraints.: */ + if (strcmp(iomodel->meshtype,"3d")==0){ + /*On a 3d mesh, we may have collapsed elements, hence dead grids. Freeze them out: */ + if (!iomodel->gridonbed[i]){ + for(k=1;k<=node_numdofs;k++){ + node->FreezeDof(k); + } + } + } + /*Add node to nodes dataset: */ + nodes->AddObject(node); + + #ifdef _PARALLEL_ + } //if((my_grids[i]==1)) + #endif + } + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + elements->Presort(); + nodes->Presort(); + materials->Presort(); + + /*Clean fetched data: */ + xfree((void**)&iomodel->deadgrids); + xfree((void**)&iomodel->x); + xfree((void**)&iomodel->y); + xfree((void**)&iomodel->z); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->gridonbed); + xfree((void**)&iomodel->gridonsurface); + xfree((void**)&iomodel->uppernodes); + xfree((void**)&iomodel->gridonicesheet); + xfree((void**)&iomodel->gridoniceshelf); + + + /*Keep partitioning information into iomodel*/ + iomodel->epart=epart; + iomodel->my_grids=my_grids; + iomodel->my_bordergrids=my_bordergrids; + + /*Free ressources:*/ + #ifdef _PARALLEL_ + xfree((void**)&all_numgrids); + xfree((void**)&npart); + VecFree(&gridborder); + #endif + + cleanup_and_return: + + /*Assign output pointer: */ + *pelements=elements; + *pnodes=nodes; + *pmaterials=materials; + +} Index: /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp (revision 2722) +++ /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateLoadsBalancedvelocities.cpp (revision 2722) @@ -0,0 +1,35 @@ +/*! \file CreateLoadsBalancedvelocities.c: + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateLoadsBalancedvelocities" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../../include/macros.h" +#include "../IoModel.h" + + +void CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel,ConstDataHandle iomodel_handle){ + + DataSet* loads = NULL; + + /*Create loads: */ + loads = new DataSet(LoadsEnum()); + + + /*All our datasets are already order by ids. Set presort flag so that later on, when sorting is requested on these + * datasets, it will not be redone: */ + loads->Presort(); + + cleanup_and_return: + + /*Assign output pointer: */ + *ploads=loads; + +} + + Index: /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp (revision 2722) +++ /issm/trunk/src/c/ModelProcessorx/Balancedvelocities/CreateParametersBalancedvelocities.cpp (revision 2722) @@ -0,0 +1,103 @@ +/*!\file: CreateParametersBalancedvelocities.cpp + * \brief driver for creating parameters dataset, for prognostic analysis. + */ + +#undef __FUNCT__ +#define __FUNCT__ "CreateParametersBalancedvelocities" + +#include "../../DataSet/DataSet.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" +#include "../../objects/objects.h" +#include "../../shared/shared.h" +#include "../IoModel.h" + +void CreateParametersBalancedvelocities(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle){ + + Param* param = NULL; + DataSet* parameters=NULL; + int count; + int i; + int dim; + + double* vx=NULL; + double* vy=NULL; + double* vz=NULL; + double* u_g=NULL; + double* pressure=NULL; + double* temperature=NULL; + double* thickness=NULL; + double* surface=NULL; + double* bed=NULL; + double* accumulation=NULL; + double* melting=NULL; + + /*recover parameters : */ + parameters=*pparameters; + + count=parameters->Size(); + + /*Get vx and vy: */ + IoModelFetchData(&vx,NULL,NULL,iomodel_handle,"vx"); + IoModelFetchData(&vy,NULL,NULL,iomodel_handle,"vy"); + IoModelFetchData(&vz,NULL,NULL,iomodel_handle,"vz"); + + u_g=(double*)xcalloc(iomodel->numberofnodes*3,sizeof(double)); + + if(vx) for(i=0;inumberofnodes;i++)u_g[3*i+0]=vx[i]/iomodel->yts; + if(vy) for(i=0;inumberofnodes;i++)u_g[3*i+1]=vy[i]/iomodel->yts; + if(vz) for(i=0;inumberofnodes;i++)u_g[3*i+2]=vz[i]/iomodel->yts; + + count++; + param= new Param(count,"u_g",DOUBLEVEC); + param->SetDoubleVec(u_g,3*iomodel->numberofnodes,3); + parameters->AddObject(param); + + + xfree((void**)&vx); + xfree((void**)&vy); + xfree((void**)&vz); + xfree((void**)&u_g); + + /*Get thickness: */ + IoModelFetchData(&thickness,NULL,NULL,iomodel_handle,"thickness"); + + count++; + param= new Param(count,"h_g",DOUBLEVEC); + if(thickness) param->SetDoubleVec(thickness,iomodel->numberofnodes,1); + else param->SetDoubleVec(thickness,0,0); + parameters->AddObject(param); + + /*Free thickness: */ + xfree((void**)&thickness); + + /*Get melting: */ + IoModelFetchData(&melting,NULL,NULL,iomodel_handle,"melting"); + if(melting) for(i=0;inumberofnodes;i++)melting[i]=melting[i]/iomodel->yts; + + count++; + param= new Param(count,"m_g",DOUBLEVEC); + if(melting) param->SetDoubleVec(melting,iomodel->numberofnodes,1); + else param->SetDoubleVec(melting,0,1); + parameters->AddObject(param); + + /*Free melting: */ + xfree((void**)&melting); + + /*Get accumulation: */ + IoModelFetchData(&accumulation,NULL,NULL,iomodel_handle,"accumulation"); + if(accumulation) for(i=0;inumberofnodes;i++)accumulation[i]=accumulation[i]/iomodel->yts; + + count++; + param= new Param(count,"a_g",DOUBLEVEC); + if(accumulation) param->SetDoubleVec(accumulation,iomodel->numberofnodes,1); + else param->SetDoubleVec(accumulation,0,0); + parameters->AddObject(param); + + /*Free accumulation: */ + xfree((void**)&accumulation); + + + /*Assign output pointer: */ + *pparameters=parameters; +} Index: /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 2721) +++ /issm/trunk/src/c/ModelProcessorx/CreateDataSets.cpp (revision 2722) @@ -94,4 +94,12 @@ } + else if (iomodel->analysis_type==BalancedvelocitiesAnalysisEnum()){ + + CreateElementsNodesAndMaterialsBalancedvelocities(pelements,pnodes,pmaterials, iomodel,iomodel_handle); + CreateConstraintsBalancedvelocities(pconstraints,iomodel,iomodel_handle); + CreateLoadsBalancedvelocities(ploads,iomodel,iomodel_handle); + CreateParametersBalancedvelocities(pparameters,iomodel,iomodel_handle); + + } else{ Index: /issm/trunk/src/c/ModelProcessorx/IoModel.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 2721) +++ /issm/trunk/src/c/ModelProcessorx/IoModel.h (revision 2722) @@ -177,23 +177,22 @@ 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; + int* penaltypartitioning; }; - /*constructor and destructor: */ IoModel* NewIoModel(void); - void DeleteIoModel( IoModel** pthis); + void DeleteIoModel( IoModel** pthis); /*Echo: */ - void IoModelEcho(IoModel* iomodel,int which_part,int rank); + void IoModelEcho(IoModel* iomodel,int which_part,int rank); /*Initialization with matlab workspace data, or marshall binary data: */ - int IoModelInit(IoModel** piomodel,ConstDataHandle iomodel_handle); + int IoModelInit(IoModel** piomodel,ConstDataHandle iomodel_handle); /*Creation of fem datasets: general drivers*/ - void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateDataSets(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, DataSet** pconstraints, DataSet** ploads,DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParameters(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); @@ -203,57 +202,62 @@ void CreateElementsNodesAndMaterialsDiagnosticHoriz(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticHoriz(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsDiagnosticHoriz(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersDiagnosticHoriz(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*diagnostic vertical*/ void CreateElementsNodesAndMaterialsDiagnosticVert(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticVert(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateLoadsDiagnosticVert(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*diagnostic hutter*/ void CreateElementsNodesAndMaterialsDiagnosticHutter(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticHutter(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateLoadsDiagnosticHutter(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*diagnostic stokes*/ void CreateElementsNodesAndMaterialsDiagnosticStokes(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsDiagnosticStokes(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateLoadsDiagnosticStokes(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*slope compute*/ void CreateElementsNodesAndMaterialsSlopeCompute(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsSlopeCompute(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateLoadsSlopeCompute(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); /*control:*/ - void CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateParametersControl(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*thermal:*/ void CreateElementsNodesAndMaterialsThermal(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsThermal(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsThermal(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersThermal(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*melting:*/ void CreateElementsNodesAndMaterialsMelting(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsMelting(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsMelting(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersMelting(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*prognostic:*/ void CreateElementsNodesAndMaterialsPrognostic(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsPrognostic(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsPrognostic(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersPrognostic(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*balancedthickness:*/ void CreateElementsNodesAndMaterialsBalancedthickness(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); void CreateConstraintsBalancedthickness(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); - void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); - void CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsBalancedthickness(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersBalancedthickness(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); + + /*balancedvelocities:*/ + void CreateElementsNodesAndMaterialsBalancedvelocities(DataSet** pelements,DataSet** pnodes, DataSet** pmaterials, IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateConstraintsBalancedvelocities(DataSet** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle); + void CreateLoadsBalancedvelocities(DataSet** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle); + void CreateParametersBalancedvelocities(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); /*qmu: */ void CreateParametersQmu(DataSet** pparameters,IoModel* iomodel,ConstDataHandle iomodel_handle); - #endif /* _IOMODEL_H */ Index: /issm/trunk/src/c/objects/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Penta.cpp (revision 2721) +++ /issm/trunk/src/c/objects/Penta.cpp (revision 2722) @@ -325,4 +325,8 @@ CreateKMatrixBalancedthickness( Kgg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==BalancedvelocitiesAnalysisEnum()){ + + CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type); + } else if (analysis_type==ThermalAnalysisEnum()){ @@ -344,4 +348,27 @@ void Penta::CreateKMatrixBalancedthickness(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreateKMatrix(Kgg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; + +} +/*}}}*/ +/*FUNCTION CreateKMatrixBalancedvelocities {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreateKMatrixBalancedvelocities" + +void Penta::CreateKMatrixBalancedvelocities(Mat Kgg,void* inputs,int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ @@ -1410,4 +1437,8 @@ CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==BalancedvelocitiesAnalysisEnum()){ + + CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type); + } else if (analysis_type==ThermalAnalysisEnum()){ @@ -1429,4 +1460,26 @@ void Penta::CreatePVectorBalancedthickness( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ + + /*Collapsed formulation: */ + Tria* tria=NULL; + + /*If on water, skip: */ + if(onwater)return; + + /*Is this element on the bed? :*/ + if(!onbed)return; + + /*Spawn Tria element from the base of the Penta: */ + tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. + tria->CreatePVector(pg,inputs, analysis_type,sub_analysis_type); + delete tria; + return; +} +/*}}}*/ +/*FUNCTION CreatePVectorBalancedvelocities {{{1*/ +#undef __FUNCT__ +#define __FUNCT__ "Penta::CreatePVectorBalancedvelocities" + +void Penta::CreatePVectorBalancedvelocities( Vec pg, void* inputs, int analysis_type,int sub_analysis_type){ /*Collapsed formulation: */ Index: /issm/trunk/src/c/objects/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Penta.h (revision 2721) +++ /issm/trunk/src/c/objects/Penta.h (revision 2722) @@ -118,4 +118,6 @@ void CreateKMatrixBalancedthickness(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); void CreatePVectorBalancedthickness( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); + void CreateKMatrixBalancedvelocities(Mat Kgg,void* vinputs,int analysis_type,int sub_analysis_type); + void CreatePVectorBalancedvelocities( Vec pg, void* vinputs, int analysis_type,int sub_analysis_type); void CreateKMatrixDiagnosticStokes( Mat Kgg, void* vinputs, int analysis_type,int sub_analysis_type); Index: /issm/trunk/src/c/objects/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Tria.cpp (revision 2721) +++ /issm/trunk/src/c/objects/Tria.cpp (revision 2722) @@ -303,4 +303,9 @@ } + else if (analysis_type==BalancedvelocitiesAnalysisEnum()){ + + CreateKMatrixBalancedvelocities( Kgg,inputs,analysis_type,sub_analysis_type); + + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet")); @@ -447,4 +452,216 @@ for( i=0; iartdiff){ + + /* Compute artificial diffusivity */ + KDL[0][0]=DL_scalar*K[0][0]; + KDL[1][1]=DL_scalar*K[1][1]; + + TripleMultiply( &Bprime[0][0],2,numdof,1, + &KDL[0][0],2,2,0, + &Bprime[0][0],2,numdof,0, + &Ke_gg_gaussian[0][0],0); + + /* Add artificial diffusivity matrix */ + for( i=0; iRecover("velocity_average",&vxvy_list[0][0],2,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find velocity_average in inputs!"); + + for(i=0;iartdiff){ + //Get the Jacobian determinant + gauss_l1l2l3[0]=1.0/3.0; gauss_l1l2l3[1]=1.0/3.0; gauss_l1l2l3[2]=1.0/3.0; + GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss_l1l2l3); + + //Build K matrix (artificial diffusivity matrix) + v_gauss[0]=1.0/3.0*(vxvy_list[0][0]+vxvy_list[1][0]+vxvy_list[2][0]); + v_gauss[1]=1.0/3.0*(vxvy_list[0][1]+vxvy_list[1][1]+vxvy_list[2][1]); + + K[0][0]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[0]); //pow should be zero!! + K[1][1]=pow(10,2)*pow(Jdettria,(double).5)/2.0*fabs(v_gauss[1]); + } + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igartdiff){ @@ -1550,4 +1767,8 @@ CreatePVectorBalancedthickness( pg,inputs,analysis_type,sub_analysis_type); } + else if (analysis_type==BalancedvelocitiesAnalysisEnum()){ + + CreatePVectorBalancedvelocities( pg,inputs,analysis_type,sub_analysis_type); + } else{ throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet")); @@ -1630,4 +1851,95 @@ /* Get accumulation, melting and thickness at gauss point */ + GetParameterValue(&accumulation_g, &accumulation_list[0],gauss_l1l2l3); + GetParameterValue(&melting_g, &melting_list[0],gauss_l1l2l3); + + /* Add value into pe_g: */ + for( i=0; iRecover("accumulation",&accumulation_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find accumulation in inputs!"); + found=inputs->Recover("melting",&melting_list[0],1,dofs,numgrids,(void**)nodes); + if(!found)throw ErrorException(__FUNCT__," could not find melting in inputs!"); + + /* Get node coordinates and dof list: */ + GetElementNodeData( &xyz_list[0][0], nodes, numgrids); + GetDofList(&doflist[0],&numberofdofspernode); + + /* Get gaussian points and weights (make this a statically initialized list of points? fstd): */ + GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); + + /* Start looping on the number of gaussian points: */ + for (ig=0; igGetFormulation(BalancedthicknessAnalysisEnum()); + } + if(analysis_type==BalancedvelocitiesAnalysisEnum()){ + fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum()); } fem_t=model->GetFormulation(ThermalAnalysisEnum()); @@ -327,4 +333,27 @@ VecFree(&h_g); } + else if(strcmp(result->GetFieldName(),"v_g")==0){ + /*easy, v_g is of size numberofnodes, on 1 dof, just repartition: */ + result->GetField(&v_g); + VecToMPISerial(&v_g_serial,v_g); + fem_p->parameters->FindParam(&numberofnodes,"numberofnodes"); + VecToMPISerial(&partition,fem_p->partition->vector); + + vel=(double*)xmalloc(numberofnodes*sizeof(double)); + + for(i=0;iSize()+1,result->GetTime(),result->GetStep(),"vel",vel,numberofnodes); + newresults->AddObject(newresult); + + /*do some cleanup: */ + xfree((void**)&v_g_serial); + xfree((void**)&vel); + xfree((void**)&partition); + VecFree(&v_g); + } else if(strcmp(result->GetFieldName(),"s_g")==0){ /*easy, s_g is of size numberofnodes, on 1 dof, just repartition: */ Index: /issm/trunk/src/c/parallel/balancedvelocities.cpp =================================================================== --- /issm/trunk/src/c/parallel/balancedvelocities.cpp (revision 2722) +++ /issm/trunk/src/c/parallel/balancedvelocities.cpp (revision 2722) @@ -0,0 +1,157 @@ +/*!\file: balancedvelocities.cpp + * \brief: balancedvelocities solution + */ + +#include "../issm.h" +#include "./parallel.h" + +#undef __FUNCT__ +#define __FUNCT__ "balancedvelocities" + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + + +int main(int argc,char* *argv){ + + /*I/O: */ + FILE* fid=NULL; + char* inputfilename=NULL; + char* outputfilename=NULL; + char* lockname=NULL; + int numberofnodes; + double waitonlock=0; + + Model* model=NULL; + + double* u_g_serial=NULL; + double* melting_g=NULL; + double* accumulation_g=NULL; + int qmu_analysis; + + /*Results: */ + DataSet* results=NULL; + DataSet* processedresults=NULL; + Result* result=NULL; + + ParameterInputs* inputs=NULL; + Param* param=NULL; + + /*time*/ + double start, finish; + double start_core, finish_core; + double start_init, finish_init; + + MODULEBOOT(); + + #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_)) + throw ErrorException(__FUNCT__," parallel executable was compiled without support of parallel libraries!"); + #endif + + /*Initialize Petsc and get start time*/ + PetscInitialize(&argc,&argv,(char *)0,""); + MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime(); + + /*Size and rank: */ + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + MPI_Comm_size(MPI_COMM_WORLD,&num_procs); + + _printf_("recover , input file name and output file name:\n"); + inputfilename=argv[2]; + outputfilename=argv[3]; + lockname=argv[4]; + + /*Initialize model structure: */ + MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime(); + model=new Model(); + + /*Open handle to data on disk: */ + fid=pfopen(inputfilename,"rb"); + + /*Initialize model structure: */ + model=new Model(); + + _printf_("read and create finite element model:\n"); + model->AddFormulation(fid,BalancedvelocitiesAnalysisEnum()); + + /*recover parameters: */ + model->FindParam(&waitonlock,"waitonlock"); + model->FindParam(&qmu_analysis,"qmu_analysis"); + + _printf_("initialize inputs:\n"); + + model->FindParam(&u_g_serial,NULL,NULL,"u_g",BalancedvelocitiesAnalysisEnum()); + model->FindParam(&melting_g,NULL,NULL,"m_g",BalancedvelocitiesAnalysisEnum()); + model->FindParam(&accumulation_g,NULL,NULL,"a_g",BalancedvelocitiesAnalysisEnum()); + model->FindParam(&numberofnodes,"numberofnodes"); + + inputs=new ParameterInputs; + inputs->Add("velocity",u_g_serial,3,numberofnodes); + inputs->Add("melting",melting_g,1,numberofnodes); + inputs->Add("accumulation",accumulation_g,1,numberofnodes); + + _printf_("initialize results:\n"); + results=new DataSet(ResultsEnum()); + MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime(); + + /*are we running the solutoin sequence, or a qmu wrapper around it? : */ + if(!qmu_analysis){ + + /*run balancedvelocities analysis: */ + _printf_("call computational core:\n"); + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + balancedvelocities_core(results,model,inputs); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + + } + else{ + + /*run qmu analysis: */ + _printf_("calling qmu analysis on balancedvelocities core:\n"); + + #ifdef _HAVE_DAKOTA_ + MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( ); + Qmux(model,inputs,BalancedvelocitiesAnalysisEnum(),NoneAnalysisEnum()); + MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( ); + #else + throw ErrorException(__FUNCT__," Dakota not present, cannot do qmu!"); + #endif + } + + /*Add analysis_type to results: */ + result=new Result(results->Size()+1,0,1,"analysis_type","balancedvelocities"); + results->AddObject(result); + + _printf_("process results:\n"); + ProcessResults(&processedresults,results,model,BalancedvelocitiesAnalysisEnum()); + + _printf_("write results to disk:\n"); + OutputResults(processedresults,outputfilename); + + if (waitonlock>0){ + _printf_("write lock file:\n"); + WriteLockFile(lockname); + } + + /*Free ressources:*/ + delete processedresults; + delete results; + delete model; + delete inputs; + + /*Get finish time and close*/ + MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( ); + _printf_("\n %-34s %f seconds \n","Model initialization elapsed time:",finish_init-start_init); + _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core); + _printf_(" %-34s %f seconds\n\n","Total elapsed time:",finish-start); + _printf_("closing MPI and Petsc\n"); + PetscFinalize(); + + /*end module: */ + MODULEEND(); + + return 0; //unix success return; +} Index: /issm/trunk/src/c/parallel/balancedvelocities_core.cpp =================================================================== --- /issm/trunk/src/c/parallel/balancedvelocities_core.cpp (revision 2722) +++ /issm/trunk/src/c/parallel/balancedvelocities_core.cpp (revision 2722) @@ -0,0 +1,64 @@ +/*!\file: balancedvelocities_core.cpp + * \brief: core of the balancedvelocities solution + */ + +#undef __FUNCT__ +#define __FUNCT__ "balancedvelocities_core" + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./parallel.h" +#include "../issm.h" + +void balancedvelocities_core(DataSet* results,Model* model,ParameterInputs* inputs){ + + extern int my_rank; + + /*output: */ + Result* result=NULL; + + /*intermediary: */ + Vec u_g=NULL; + + /*solutions: */ + Vec v_g=NULL; + + /*flags: */ + int verbose=0; + int numberofdofspernode; + int numberofnodes; + int dofs[2]={1,1}; + + /*fem balancedvelocities model: */ + FemModel* fem_p=NULL; + + + /*recover fem model: */ + fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum()); + + //first recover parameters common to all solutions + model->FindParam(&verbose,"verbose"); + model->FindParam(&numberofnodes,"numberofnodes"); + model->FindParam(&numberofdofspernode,"numberofdofspernode"); + + _printf_("depth averaging velocity...\n"); + u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity + FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"velocity"); + inputs->Add("velocity_average",u_g,2,numberofnodes); + + _printf_("call computational core:\n"); + diagnostic_core_linear(&v_g,fem_p,inputs,BalancedvelocitiesAnalysisEnum(),NoneAnalysisEnum()); + + _printf_("extrude computed thickness on all layers:\n"); + FieldExtrudex( v_g, fem_p->elements,fem_p->nodes, fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0); + + /*Plug results into output dataset: */ + result=new Result(results->Size()+1,0,1,"v_g",v_g); + results->AddObject(result); + + /*Free ressources:*/ + VecFree(&u_g); + VecFree(&v_g); +} Index: /issm/trunk/src/c/parallel/parallel.h =================================================================== --- /issm/trunk/src/c/parallel/parallel.h (revision 2721) +++ /issm/trunk/src/c/parallel/parallel.h (revision 2722) @@ -18,4 +18,5 @@ void prognostic_core(DataSet* results,Model* model, ParameterInputs* inputs); void balancedthickness_core(DataSet* results,Model* model, ParameterInputs* inputs); +void balancedvelocities_core(DataSet* results,Model* model, ParameterInputs* inputs); void control_core(DataSet* results,Model* model, ParameterInputs* inputs); Index: /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 2721) +++ /issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp (revision 2722) @@ -59,4 +59,7 @@ numdofs=1; } + else if (analysis_type==BalancedvelocitiesAnalysisEnum()){ + numdofs=1; + } else throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis type: ",analysis_type," not implemented yet!")); Index: /issm/trunk/src/m/classes/public/process_solve_options.m =================================================================== --- /issm/trunk/src/m/classes/public/process_solve_options.m (revision 2721) +++ /issm/trunk/src/m/classes/public/process_solve_options.m (revision 2722) @@ -17,9 +17,10 @@ %check solution type is supported -if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient','balancedthickness'}), +if ~ismemberi(analysis_type,{'diagnostic','prognostic','thermal','steadystate','parameters','transient',... + 'balancedthickness','balancedvelocities'}), error(['process_solve_options error message: analysis_type ' analysis_type ' not supported yet!']); else %convert to enum - outoptions.analysis_type=eval([upper(analysis_type(1)) lower(analysis_type(2:end)) 'AnalysisEnum']); + outoptions.analysis_type=eval([upper(analysis_type(1)) lower(analysis_type(2:end)) 'AnalysisEnum()']); end if ~ismemberi(sub_analysis_type,{'steady','transient','none','horiz','adjoint','gradient','inverse','vert',''}), @@ -27,5 +28,5 @@ else %convert to enum - outoptions.sub_analysis_type=eval([upper(sub_analysis_type(1)) lower(sub_analysis_type(2:end)) 'AnalysisEnum']); + outoptions.sub_analysis_type=eval([upper(sub_analysis_type(1)) lower(sub_analysis_type(2:end)) 'AnalysisEnum()']); end Index: /issm/trunk/src/m/classes/public/solve.m =================================================================== --- /issm/trunk/src/m/classes/public/solve.m (revision 2721) +++ /issm/trunk/src/m/classes/public/solve.m (revision 2722) @@ -73,4 +73,8 @@ elseif md.analysis_type==BalancedthicknessAnalysisEnum, md=balancedthickness(md); + +elseif md.analysis_type==BalancedvelocitiesAnalysisEnum, + md=balancedvelocities(md); + else error('solution type not supported for this model configuration.'); Index: /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m =================================================================== --- /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m (revision 2721) +++ /issm/trunk/src/m/enum/AnalysisTypeFromEnum.m (revision 2722) @@ -89,4 +89,8 @@ end +if enum==BalancedvelocitiesAnalysisEnum(), + string='balancedvelocities'; +end + if enum==MeltingAnalysisEnum(), string='melting'; Index: /issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m =================================================================== --- /issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m (revision 2722) +++ /issm/trunk/src/m/enum/BalancedvelocitiesAnalysisEnum.m (revision 2722) @@ -0,0 +1,9 @@ +function macro=BalancedvelocitiesAnalysisEnum() +%BALANCEDVELOCITIESANALYSISENUM - Enum of BalancedvelocitiesAnalysis +% +% file generated by src/c/EnumDefinitions/SynchronizeMatlabEnum +% +% Usage: +% macro=BalancedvelocitiesAnalysisEnum() + +macro=252; Index: /issm/trunk/src/m/solutions/jpl/balancedvelocities.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedvelocities.m (revision 2722) +++ /issm/trunk/src/m/solutions/jpl/balancedvelocities.m (revision 2722) @@ -0,0 +1,37 @@ +function md=balancedvelocities(md) +%BALANCEDVELOCITIES - balancedvelocities solution sequence. +% +% Usage: +% md=balancedvelocities(md) + + %timing + t1=clock; + + %Build all models requested for diagnostic simulation + models.analysis_type=BalancedvelocitiesAnalysisEnum; %needed for processresults + + displaystring(md.verbose,'%s',['reading balancedvelocities model data']); + md.analysis_type=BalancedvelocitiesAnalysisEnum; models.p=CreateFemModel(md); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + displaystring(md.verbose,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes); + inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double'); + + displaystring(md.verbose,'\n%s',['call computational core:']); + results=balancedvelocities_core(models,inputs,BalancedvelocitiesAnalysisEnum(),NoneAnalysisEnum()); + + displaystring(md.verbose,'\n%s',['load results...']); + if ~isstruct(md.results), md.results=struct(); end + md.results.balancedvelocities=processresults(models, results); + + %stop timing + t2=clock; + displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m (revision 2722) +++ /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m (revision 2722) @@ -0,0 +1,24 @@ +function results=balancedvelocities_core(models,inputs,analysis_type,sub_analysis_type) +%BALANCEDVELOCITIES_CORE - linear solution sequence +% +% Usage: +% v_g=balancedvelocities_core(m,inputs,analysis_type,sub_analysis_type) + + %get FE model + m=models.p; + results.time=0; + results.step=1; + + displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']); + %Take only the first two dofs of m.parameters.u_g + u_g=get(inputs,'velocity',[1 1 0 0]); + u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity'); + inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes); + + displaystring(m.parameters.verbose,'\n%s',['call computational core:']); + results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type); + + displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']); + results.v_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0); + +end %end function Index: /issm/trunk/src/m/solutions/jpl/processresults.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/processresults.m (revision 2721) +++ /issm/trunk/src/m/solutions/jpl/processresults.m (revision 2722) @@ -126,4 +126,11 @@ newresults(i).parameter=param_g; + elseif strcmpi(resultsname{j},'v_g'), + + v_g=results(i).v_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).vel=v_g; + else