Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 652) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 653) @@ -31,4 +31,5 @@ double* dirichletvalues_diag=NULL; double* gridondirichlet_diag=NULL; + double* gridonhutter=NULL; /*Create constraints: */ @@ -41,4 +42,5 @@ ModelFetchData((void**)&gridondirichlet_diag,NULL,NULL,model_handle,"gridondirichlet_diag","Matrix","Mat"); ModelFetchData((void**)&dirichletvalues_diag,NULL,NULL,model_handle,"dirichletvalues_diag","Matrix","Mat"); + ModelFetchData((void**)&gridonhutter,NULL,NULL,model_handle,"gridonhutter","Matrix","Mat"); count=0; @@ -51,5 +53,5 @@ #endif - if ((int)gridondirichlet_diag[i]){ + if ((int)gridondirichlet_diag[i] | (int)gridonhutter[i]){ /*This grid needs to be spc'd to vx_obs and vy_obs:*/ @@ -87,4 +89,5 @@ xfree((void**)&gridondirichlet_diag); xfree((void**)&dirichletvalues_diag); + xfree((void**)&gridonhutter); cleanup_and_return: Index: /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 652) +++ /issm/trunk/src/c/ModelProcessorx/DiagnosticHoriz/CreateElementsNodesAndMaterialsDiagnosticHoriz.cpp (revision 653) @@ -229,4 +229,5 @@ 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->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"); @@ -241,89 +242,91 @@ #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]; - - /*meling and accumulation*/ - tria_melting[0]=*(model->melting+ ((int)*(model->elements+elements_width*i+0)-1)); - tria_melting[1]=*(model->melting+ ((int)*(model->elements+elements_width*i+1)-1)); - tria_melting[2]=*(model->melting+ ((int)*(model->elements+elements_width*i+2)-1)); - - tria_accumulation[0]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+0)-1)); - tria_accumulation[1]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+1)-1)); - tria_accumulation[2]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+2)-1)); - - /*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_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff); - - /*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 + if (*(model->elements_type+2*i+0)==MacAyealEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements. + + /*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]; + + /*meling and accumulation*/ + tria_melting[0]=*(model->melting+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_melting[1]=*(model->melting+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_melting[2]=*(model->melting+ ((int)*(model->elements+elements_width*i+2)-1)); + + tria_accumulation[0]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+0)-1)); + tria_accumulation[1]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+1)-1)); + tria_accumulation[2]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+2)-1)); + + /*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_melting,tria_accumulation,tria_geothermalflux,tria_friction_type, tria_p, tria_q, tria_shelf, tria_meanvel, tria_epsvel, tria_viscosity_overshoot,tria_artdiff); + + /*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 + } //if(!HutterEnum) #ifdef _PARALLEL_ @@ -345,4 +348,5 @@ xfree((void**)&model->B); xfree((void**)&model->n); + xfree((void**)&model->elements_type); } @@ -372,84 +376,86 @@ #endif + if (*(model->elements_type+2*i+0)==MacAyealEnum() | *(model->elements_type+2*i+0)==PattynEnum()){ //elements of type 1 are Hutter type Tria. Don't create this elements. - /*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)); - penta_melting[j]=*(model->melting+ ((int)*(model->elements+elements_width*i+j)-1)); - penta_accumulation[j]=*(model->accumulation+ ((int)*(model->elements+elements_width*i+j)-1)); + /*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)); + penta_melting[j]=*(model->melting+ ((int)*(model->elements+elements_width*i+j)-1)); + penta_accumulation[j]=*(model->accumulation+ ((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,penta_stokesreconditioning); + + /*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 } - - /*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,penta_stokesreconditioning); - - /*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_ @@ -572,5 +578,4 @@ ModelFetchData((void**)&model->gridonsurface,NULL,NULL,model_handle,"gridonsurface","Matrix","Mat"); - /*Get number of dofs per node: */ DistributeNumDofs(&node_numdofs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/c/ModelProcessorx/Model.cpp =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 652) +++ /issm/trunk/src/c/ModelProcessorx/Model.cpp (revision 653) @@ -53,4 +53,6 @@ model->uppernodes=NULL; model->gridonhutter=NULL; + model->gridonmacayeal=NULL; + model->gridonpattyn=NULL; model->vx_obs=NULL; @@ -203,8 +205,10 @@ xfree((void**)&model->elements_type); xfree((void**)&model->gridonhutter); + xfree((void**)&model->gridonmacayeal); if (strcmp(model->meshtype,"3d")==0){ xfree((void**)&model->elements2d); xfree((void**)&model->deadgrids); xfree((void**)&model->uppernodes); + xfree((void**)&model->gridonpattyn); } xfree((void**)&model->solverstring); Index: /issm/trunk/src/c/ModelProcessorx/Model.h =================================================================== --- /issm/trunk/src/c/ModelProcessorx/Model.h (revision 652) +++ /issm/trunk/src/c/ModelProcessorx/Model.h (revision 653) @@ -42,4 +42,6 @@ int isstokes; double* gridonhutter; + double* gridonmacayeal; + double* gridonpattyn; /*results: */