Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.cpp (revision 15484) @@ -6,22 +6,24 @@ #include "../../shared/shared.h" #include "../../toolkits/toolkits.h" +#include "../ModelProcessorx/ModelProcessorx.h" -void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type){ +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element){ /*intermediary: */ - int i,j; - IssmDouble yts; - bool transient = false; + int i,j,count; + IssmDouble yts; FILE *fid = NULL; - int code = 0; - int vector_layout = 0; + int code,vector_layout; IssmDouble *times = NULL; IssmDouble *values = NULL; bool spcpresent = false; - int count = 0; + + /*P2 finite elements*/ + int v1,v2; + bool *my_edges = NULL; /*variables being fetched: */ - IssmDouble *IssmDoublevector = NULL; - int M,N; + IssmDouble *spcdata = NULL; + int M,N; /*Fetch parameters: */ @@ -30,26 +32,36 @@ /*First of, find the record for the enum, and get code of data type: */ fid=iomodel->SetFilePointerToData(&code, &vector_layout,vector_enum); - if(code!=7)_error_("expecting a IssmDouble vector for constraints with enum " << EnumToStringx(vector_enum)); if(vector_layout!=1)_error_("expecting a nodal vector for constraints with enum " << EnumToStringx(vector_enum)); /*Fetch vector:*/ - iomodel->FetchData(&IssmDoublevector,&M,&N,vector_enum); + iomodel->FetchData(&spcdata,&M,&N,vector_enum); - /*Transient or static?:*/ + /*Partition edges if we are using P2 finite elements*/ + if(finite_element==P2Enum){ + EdgesPartitioning(&my_edges,iomodel); + } + if(M==iomodel->numberofvertices){ - /*static: just create Constraints objects*/ + /*Static constraint*/ count=0; - - /*Create Constraints from x,y,z: */ for (i=0;inumberofvertices;i++){ - - /*keep only this partition's nodes:*/ if((iomodel->my_vertices[i])){ - - if (!xIsNan(IssmDoublevector[i])){ - - constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,IssmDoublevector[i],analysis_type)); + if (!xIsNan(spcdata[i])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,spcdata[i],analysis_type)); count++; + } + } + } + if(finite_element==P2Enum){ + for(i=0;inumberofedges;i++){ + if(my_edges[i]){ + v1 = iomodel->edges[4*i+0]-1; + v2 = iomodel->edges[4*i+1]-1; + if(!xIsNan(spcdata[v1]) && !xIsNan(spcdata[v2])){ + constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, + 1,(spcdata[v1]+spcdata[v2])/2.,analysis_type)); + count++; + } } } @@ -63,12 +75,8 @@ /*figure out times: */ times=xNew(N); - for(j=0;jnumberofvertices;i++){ - - /*keep only this partition's nodes:*/ if((iomodel->my_vertices[i])){ @@ -77,5 +85,5 @@ spcpresent=false; for(j=0;j(values[j]))spcpresent=true; //NaN means no spc by default } @@ -88,4 +96,23 @@ } } + if(finite_element==P2Enum){ + for(i=0;inumberofedges;i++){ + if(my_edges[i]){ + v1 = iomodel->edges[4*i+0]-1; + v2 = iomodel->edges[4*i+1]-1; + values=xNew(N); + spcpresent=false; + for(j=0;j(values[j])) spcpresent=true; //NaN means no spc by default + } + if(spcpresent){ + constraints->AddObject(new SpcTransient(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,1, + N,times,values,analysis_type)); + count++; + } + } + } + } } else{ @@ -94,5 +121,5 @@ /*Free ressources:*/ - xDelete(IssmDoublevector); + xDelete(spcdata); xDelete(times); xDelete(values); Index: /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15483) +++ /issm/trunk-jpl/src/c/modules/IoModelToConstraintsx/IoModelToConstraintsx.h (revision 15484) @@ -8,5 +8,5 @@ /* local prototypes: */ -void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type); +void IoModelToConstraintsx(Constraints* constraints,IoModel* iomodel,int vector_enum,int analysis_type,int finite_element); #endif /* _IOMODELTOELEMENTINPUTX_H */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp (revision 15484) @@ -16,5 +16,5 @@ /*Do not add constraints in DG*/ if(stabilization!=3){ - IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,BalancethicknessSpcthicknessEnum,BalancethicknessAnalysisEnum,P1Enum); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15484) @@ -11,16 +11,16 @@ /*Intermediary*/ - int i,j; - int count; - IssmDouble yts; - IssmDouble g; - IssmDouble rho_ice; - IssmDouble stokesreconditioning; - bool isstokes,isl1l2,ismacayealpattyn; - int fe_ssa; - bool spcpresent=false; - int Mx,Nx; - int My,Ny; - int Mz,Nz; + int i,j; + int count; + IssmDouble yts; + IssmDouble g; + IssmDouble rho_ice; + IssmDouble stokesreconditioning; + bool isstokes,isl1l2,ismacayealpattyn; + int fe_ssa; + bool spcpresent = false; + int Mx,Nx; + int My,Ny; + int Mz,Nz; IssmDouble *spcvx = NULL; IssmDouble *spcvy = NULL; @@ -329,10 +329,10 @@ if(!xIsNan(spcvx[v1]) && !xIsNan(spcvx[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 1,(spcvx[v1]+spcvx[v2])/yts,DiagnosticHorizAnalysisEnum)); + 1,(spcvx[v1]+spcvx[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); count++; } if(!xIsNan(spcvy[v1]) && !xIsNan(spcvy[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 2,(spcvy[v1]+spcvy[v2])/yts,DiagnosticHorizAnalysisEnum)); + 2,(spcvy[v1]+spcvy[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); count++; } @@ -340,5 +340,5 @@ if(!xIsNan(spcvz[v1]) && !xIsNan(spcvz[v2])){ constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1, - 3,(spcvz[v1]+spcvz[v2])/yts,DiagnosticHorizAnalysisEnum)); + 3,(spcvz[v1]+spcvz[v2])/(2.*yts),DiagnosticHorizAnalysisEnum)); count++; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCEfficient/CreateConstraintsHydrologyDCEfficient.cpp (revision 15484) @@ -24,5 +24,5 @@ if(!isefficientlayer) return; - IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpceplHeadEnum,HydrologyDCEfficientAnalysisEnum,P1Enum); /*Assign output pointer: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyDCInefficient/CreateConstraintsHydrologyDCInefficient.cpp (revision 15484) @@ -19,5 +19,5 @@ if(hydrology_model!=HydrologydcEnum) return; - IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,HydrologydcSpcsedimentHeadEnum,HydrologyDCInefficientAnalysisEnum,P1Enum); /*Assign output pointer: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/HydrologyShreve/CreateConstraintsHydrologyShreve.cpp (revision 15484) @@ -20,5 +20,5 @@ if(hydrology_model!=HydrologyshreveEnum) return; - IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,HydrologyshreveSpcwatercolumnEnum,HydrologyShreveAnalysisEnum,P1Enum); /*Assign output pointer: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp (revision 15484) @@ -16,5 +16,5 @@ /*Do not add constraints in DG, they are weakly imposed*/ if(stabilization!=3){ - IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,PrognosticSpcthicknessEnum,PrognosticAnalysisEnum,P1Enum); } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 15483) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp (revision 15484) @@ -16,5 +16,5 @@ /*Only 3d mesh supported*/ if(iomodel->dim==3){ - IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum); + IoModelToConstraintsx(constraints,iomodel,ThermalSpctemperatureEnum,ThermalAnalysisEnum,P1Enum); }