Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 13050) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 13051) @@ -203,4 +203,5 @@ ./shared/Numerics/isnan.h\ ./shared/Numerics/isnan.cpp\ + ./shared/Numerics/cubic.cpp\ ./shared/Numerics/extrema.cpp\ ./shared/Numerics/XZvectorsToCoordinateSystem.cpp\ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.cpp (revision 13051) @@ -1677,30 +1677,33 @@ //Need to know the type of approximation for this element if(iomodel->Data(FlowequationElementEquationEnum)){ - if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealApproximationEnum){ + if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealPattynApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealPattynApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==HutterApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HutterApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==StokesApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){ + this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum)); + } + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==StokesApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealStokesApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealStokesApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynStokesApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynStokesApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum)); } - else if (*(iomodel->Data(FlowequationElementEquationEnum)+index)==NoneApproximationEnum){ + else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); } else{ - _error_("Approximation type " << EnumToStringx((int)*(iomodel->Data(FlowequationElementEquationEnum)+index)) << " not supported yet"); + _error_("Approximation type " << EnumToStringx((int)iomodel->Data(FlowequationElementEquationEnum)[index]) << " not supported yet"); } } @@ -2801,14 +2804,14 @@ /*Intermediaries*/ - IssmInt i,j; - int penta_type; - int penta_node_ids[6]; - int penta_vertex_ids[6]; - IssmDouble nodeinputs[6]; - IssmDouble yts; - int stabilization; - bool dakota_analysis; - bool isstokes; - IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat; + IssmInt i,j; + int penta_type; + int penta_node_ids[6]; + int penta_vertex_ids[6]; + IssmDouble nodeinputs[6]; + IssmDouble yts; + int stabilization; + bool dakota_analysis; + bool isstokes; + IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat; /*Fetch parameters: */ @@ -4232,9 +4235,9 @@ const int numdof=NDOF1*NUMVERTICES; - int i; - int* doflist=NULL; - IssmDouble values[numdof]; - IssmDouble temp; - GaussPenta *gauss=NULL; + int i; + int *doflist = NULL; + IssmDouble values[numdof]; + IssmDouble temp; + GaussPenta *gauss = NULL; /*Get dof list: */ @@ -6200,4 +6203,6 @@ case MacAyealApproximationEnum: return CreateKMatrixDiagnosticMacAyeal2d(); + case L1L2ApproximationEnum: + return CreateKMatrixDiagnosticL1L2(); case PattynApproximationEnum: return CreateKMatrixDiagnosticPattyn(); @@ -6226,7 +6231,7 @@ /*Intermediaries*/ - int connectivity[2]; - int i,i0,i1,j0,j1; - IssmDouble one0,one1; + int connectivity[2]; + int i,i0,i1,j0,j1; + IssmDouble one0,one1; /*Initialize Element matrix*/ @@ -6453,4 +6458,103 @@ delete Ke2; delete Ke3; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticL1L2{{{*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixDiagnosticL1L2Viscous(); + ElementMatrix* Ke2=CreateKMatrixDiagnosticL1L2Friction(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticL1L2Viscous{{{*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2Viscous(void){ + + /*Constants*/ + const int numdof2d=2*NUMVERTICES2D; + + /*Intermediaries */ + int i,j; + IssmDouble Jdet,viscosity; + IssmDouble epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/ + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble B[3][numdof2d]; + IssmDouble Bprime[3][numdof2d]; + IssmDouble Ke_gg_gaussian[numdof2d][numdof2d]; + IssmDouble D[3][3]= {0.0}; // material matrix, simple scalar matrix. + Tria *tria = NULL; + Penta *pentabase = NULL; + GaussPenta *gauss = NULL; + GaussTria *gauss_tria = NULL; + + /*Find penta on bed as this is a macayeal elements: */ + pentabase=GetBasalElement(); + tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + + /*Initialize Element matrix*/ + ElementMatrix* Ke=new ElementMatrix(tria->nodes,NUMVERTICES2D,this->parameters,L1L2ApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* surf_input=inputs->GetInput(SurfaceEnum); _assert_(surf_input); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(5,5); + gauss_tria=new GaussTria(); + for(int ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + gauss->SynchronizeGaussTria(gauss_tria); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + tria->GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss_tria); + tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria); + + /*Get viscosity for L1L2 model*/ + GetL1L2Viscosity(&viscosity,&xyz_list[0][0],gauss,vx_input,vy_input,surf_input); + + for(i=0;i<3;i++) D[i][i]=2*viscosity*gauss->weight*Jdet; + + TripleMultiply( &B[0][0],3,numdof2d,1, + &D[0][0],3,3,0, + &Bprime[0][0],3,numdof2d,0, + &Ke_gg_gaussian[0][0],0); + for(i=0;ivalues[i*numdof2d+j]+=Ke_gg_gaussian[i][j]; + } + + /*Transform Coordinate System*/ + TransformStiffnessMatrixCoord(Ke,tria->nodes,NUMVERTICES2D,XYEnum); + + /*Clean up and return*/ + delete tria->matice; + delete tria; + delete gauss_tria; + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixDiagnosticL1L2Friction{{{*/ +ElementMatrix* Penta::CreateKMatrixDiagnosticL1L2Friction(void){ + + /*Initialize Element matrix and return if necessary*/ + if(IsFloating() || !IsOnBed()) return NULL; + + /*Build a tria element using the 3 nodes of the base of the penta. Then use + * the tria functionality to build a friction stiffness matrix on these 3 + * nodes: */ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction(); + delete tria->matice; delete tria; + + /*clean-up and return*/ return Ke; } @@ -7193,4 +7297,6 @@ case PattynApproximationEnum: return CreatePVectorDiagnosticPattyn(); + case L1L2ApproximationEnum: + return CreatePVectorDiagnosticL1L2(); case HutterApproximationEnum: return NULL; @@ -7346,4 +7452,18 @@ /*FUNCTION Penta::CreatePVectorDiagnosticMacAyeal{{{*/ ElementVector* Penta::CreatePVectorDiagnosticMacAyeal(void){ + + if (!IsOnBed()) return NULL; + + /*Call Tria function*/ + Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria. + ElementVector* pe=tria->CreatePVectorDiagnosticMacAyeal(); + delete tria->matice; delete tria; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorDiagnosticL1L2{{{*/ +ElementVector* Penta::CreatePVectorDiagnosticL1L2(void){ if (!IsOnBed()) return NULL; @@ -7892,10 +8012,10 @@ const int numdof=NDOF2*NUMVERTICES; - int i; - int approximation; - int* doflist=NULL; - IssmDouble vx,vy; - IssmDouble values[numdof]; - GaussPenta* gauss; + int i; + int approximation; + int *doflist = NULL; + IssmDouble vx,vy; + IssmDouble values[numdof]; + GaussPenta *gauss; /*Get approximation enum and dof list: */ @@ -8043,4 +8163,88 @@ } /*}}}*/ +/*FUNCTION Penta::GetL1L2Viscosity{{{*/ +void Penta::GetL1L2Viscosity(IssmDouble* pviscosity,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input,Input* surface_input){ + /*Compute the L1L2 viscosity + * + * 1 + * mu = - A^-1 (sigma'_e)^(1-n) + * 2 + * + * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18) + * + * L1L2 assumptions: + * + * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//| + * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2 + * + * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|) + * + * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0 */ + + _error_("Not supported yet"); + int i; + IssmDouble a,c,d,z,s,viscosity; + IssmDouble tau_perp,tau_par,eps_b,A; + IssmDouble epsilonvx[5]; /*exx eyy exy exz eyz*/ + IssmDouble epsilonvy[5]; /*exx eyy exy exz eyz*/ + IssmDouble epsilon[5]; /*exx eyy exy exz eyz*/ + IssmDouble z_list[NUMVERTICES]; + IssmDouble slope[3]; + + /*Check that both inputs have been found*/ + if (!vx_input || !vy_input || !surface_input) _error_("Input missing"); + + /*Get tau_perp*/ + for(i=0;iGetInputValue(&s,gauss); + surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss); + PentaRef::GetInputValue(&z,&z_list[0],gauss); + tau_perp = matpar->GetRhoIce() * matpar->GetG() * (s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]); + + /* Get eps_b*/ + vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss); + vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss); + for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i]; + eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]); + if(eps_b==0.){ + *pviscosity=2.5*pow(10.,17.); + return; + } + + /*Get A*/ + _assert_(matice->GetN()==3.0); + A=matice->GetA(); + + /*Solve for tau_par (http://en.wikipedia.org/wiki/Cubic_function)*/ + a=A; + c=A*tau_perp*tau_perp; + d=-eps_b; + tau_par = -1./(3.*a) * pow(0.5*(27.*a*a*d + sqrt(27.*a*a*d*27.*a*a*d -4.*pow(-3.*a*c,3.))),1./3.); + + //printf("=======================================\n"); + //printf("tau_par = %g (%g=0?)\n",tau_par,a*tau_par*tau_par*tau_par+c*tau_par+d); + //printf("a=%g c=%g d=%g\n",a,c,d); + //IssmDouble coeff[4]; + //coeff[0]=d; + //coeff[1]=c; + //coeff[2]=0.; + //coeff[3]=a; + //int numroots; + //IssmDouble roots[3]; + //cubic(coeff,roots,&numroots); + //tau_par=roots[0]; + //for(i=0;i(doflistm); xDelete(doflists); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticL1L2 {{{*/ +void Penta::InputUpdateFromSolutionDiagnosticL1L2(IssmDouble* solution){ + + const int numdof=NDOF2*NUMVERTICES; + + int i; + IssmDouble rho_ice,g; + IssmDouble values[numdof]; + IssmDouble vx[NUMVERTICES]; + IssmDouble vy[NUMVERTICES]; + IssmDouble vz[NUMVERTICES]; + IssmDouble vel[NUMVERTICES]; + IssmDouble pressure[NUMVERTICES]; + IssmDouble surface[NUMVERTICES]; + IssmDouble xyz_list[NUMVERTICES][3]; + int *doflist = NULL; + Penta *penta = NULL; + + /*Get dof list: */ + GetDofList(&doflist,L1L2ApproximationEnum,GsetEnum); + + /*Use the dof list to index into the solution vector: */ + for(i=0;i(vx[i])) _error_("NaN found in solution vector"); + if(xIsNan(vy[i])) _error_("NaN found in solution vector"); + } + + /*Get parameters fro pressure computation*/ + rho_ice=matpar->GetRhoIce(); + g=matpar->GetG(); + + /*Start looping over all elements above current element and update all inputs*/ + penta=this; + for(;;){ + + /*Get node data: */ + GetVerticesCoordinates(&xyz_list[0][0],penta->nodes,NUMVERTICES); + + /*Now Compute vel*/ + GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0 + for(i=0;iinputs->ChangeEnum(VxEnum,VxPicardEnum); + penta->inputs->ChangeEnum(VyEnum,VyPicardEnum); + penta->inputs->ChangeEnum(PressureEnum,PressurePicardEnum); + + /*Add vx and vy as inputs to the tria element: */ + penta->inputs->AddInput(new PentaP1Input(VxEnum,vx)); + penta->inputs->AddInput(new PentaP1Input(VyEnum,vy)); + penta->inputs->AddInput(new PentaP1Input(VelEnum,vel)); + penta->inputs->AddInput(new PentaP1Input(PressureEnum,pressure)); + + /*Stop if we have reached the surface*/ + if (penta->IsOnSurface()) break; + + /* get upper Penta*/ + penta=penta->GetUpperElement(); _assert_(penta->Id()!=this->id); + } + + /*Free ressources:*/ + xDelete(doflist); } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 13050) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Penta.h (revision 13051) @@ -227,4 +227,7 @@ ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void); ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void); + ElementMatrix* CreateKMatrixDiagnosticL1L2(void); + ElementMatrix* CreateKMatrixDiagnosticL1L2Viscous(void); + ElementMatrix* CreateKMatrixDiagnosticL1L2Friction(void); ElementMatrix* CreateKMatrixDiagnosticPattyn(void); ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void); @@ -245,4 +248,5 @@ void InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong); void InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong); + void InputUpdateFromSolutionDiagnosticL1L2( IssmDouble* solutiong); void InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong); void InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong); @@ -265,4 +269,5 @@ ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void); ElementVector* CreatePVectorDiagnosticMacAyealStokes(void); + ElementVector* CreatePVectorDiagnosticL1L2(void); ElementVector* CreatePVectorDiagnosticPattyn(void); ElementVector* CreatePVectorDiagnosticPattynStokes(void); @@ -273,4 +278,5 @@ ElementVector* CreatePVectorDiagnosticVertVolume(void); ElementVector* CreatePVectorDiagnosticVertBase(void); + void GetL1L2Viscosity(IssmDouble*, IssmDouble*, GaussPenta*, Input*, Input*, Input*); #endif Index: /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/classes/objects/Elements/Tria.cpp (revision 13051) @@ -232,8 +232,8 @@ /*Intermediaries */ int i,j,ig; - IssmDouble heatcapacity,latentheat; - IssmDouble Jdet,D_scalar; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble L[3]; + IssmDouble heatcapacity,latentheat; + IssmDouble Jdet,D_scalar; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble L[3]; GaussTria *gauss=NULL; @@ -1444,13 +1444,13 @@ /*Intermediaries*/ - int i,j; - int tria_vertex_ids[3]; + int i,j; + int tria_vertex_ids[3]; IssmDouble nodeinputs[3]; IssmDouble cmmininputs[3]; IssmDouble cmmaxinputs[3]; - bool control_analysis=false; - int num_control_type; + bool control_analysis = false; + int num_control_type; IssmDouble yts; - int num_cm_responses; + int num_cm_responses; /*Get parameters: */ Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.cpp (revision 13051) @@ -128,4 +128,18 @@ } /*}}}*/ +/*FUNCTION Matice::GetA {{{*/ +IssmDouble Matice::GetA(){ + /* + * A = 1/B^n + */ + + IssmDouble B,n; + + inputs->GetInputAverage(&B,MaterialsRheologyBEnum); + n=this->GetN(); + + return pow(B,-n); +} +/*}}}*/ /*FUNCTION Matice::GetB {{{*/ IssmDouble Matice::GetB(){ @@ -236,5 +250,5 @@ if(A==0){ /*Maxiviscositym viscosity for 0 shear areas: */ - viscosity=2.5*pow((IssmDouble)10,(IssmDouble)17); + viscosity=2.5*pow(10.,17.); } else{ Index: /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h (revision 13050) +++ /issm/trunk-jpl/src/c/classes/objects/Materials/Matice.h (revision 13051) @@ -58,15 +58,16 @@ /*}}}*/ /*Matice Numerics: {{{*/ - void SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin); - void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); - void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); - void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); - void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); - void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); - void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); + void GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon); + void GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon); + void GetViscosity3dStokes(IssmDouble* pviscosity3d, IssmDouble* epsilon); + void GetViscosityComplement(IssmDouble* pviscosity_complement, IssmDouble* pepsilon); + void GetViscosityDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + void GetViscosity2dDerivativeEpsSquare(IssmDouble* pmu_prime, IssmDouble* pepsilon); + IssmDouble GetA(); IssmDouble GetB(); IssmDouble GetBbar(); IssmDouble GetN(); - bool IsInput(int name); + bool IsInput(int name); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/objects/Node.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/objects/Node.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/classes/objects/Node.cpp (revision 13051) @@ -92,4 +92,7 @@ for(k=1;k<=gsize;k++) this->FreezeDof(k); } + if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==L1L2ApproximationEnum && !iomodel->Data(MeshVertexonbedEnum)[io_index]){ + for(k=1;k<=gsize;k++) this->FreezeDof(k); + } if (iomodel->Data(FlowequationVertexEquationEnum)[io_index]==MacAyealPattynApproximationEnum && iomodel->Data(FlowequationBordermacayealEnum)[io_index]){ if(!iomodel->Data(MeshVertexonbedEnum)[io_index]){ @@ -114,5 +117,5 @@ /*Diagnostic Hutter*/ if (analysis_type==DiagnosticHutterAnalysisEnum){ - if (!iomodel->Data(FlowequationVertexEquationEnum)) _error_("iomodel->vertices_type is NULL"); + _assert_(iomodel->Data(FlowequationVertexEquationEnum)); /*Constrain all nodes that are not Hutter*/ if (reCast(iomodel->Data(FlowequationVertexEquationEnum)[io_index])!=HutterApproximationEnum){ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 13051) @@ -36,4 +36,5 @@ parameters->AddObject(iomodel->CopyConstantObject(FlowequationIshutterEnum)); parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsmacayealpattynEnum)); + parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsl1l2Enum)); parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsstokesEnum)); parameters->AddObject(iomodel->CopyConstantObject(SettingsOutputFrequencyEnum)); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 13051) @@ -20,5 +20,5 @@ double rho_ice; double stokesreconditioning; - bool isstokes,ismacayealpattyn; + bool isstokes,isl1l2,ismacayealpattyn; bool spcpresent=false; int Mx,Nx; @@ -56,4 +56,5 @@ iomodel->Constant(&stokesreconditioning,DiagnosticStokesreconditioningEnum); iomodel->Constant(&isstokes,FlowequationIsstokesEnum); + iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); @@ -65,5 +66,5 @@ /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if (!ismacayealpattyn & !isstokes){ + if(!ismacayealpattyn & !isstokes & !isl1l2){ *pconstraints=constraints; return; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateLoadsDiagnosticHoriz.cpp (revision 13051) @@ -94,4 +94,8 @@ count++; } + else if ((int)*(elements_type+element)==(L1L2ApproximationEnum)){ + loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,L1L2IceFrontEnum,DiagnosticHorizAnalysisEnum)); + count++; + } else if ((int)*(elements_type+element)==(StokesApproximationEnum)){ loads->AddObject(new Icefront(iomodel->loadcounter+count+1,i,iomodel,StokesIceFrontEnum,DiagnosticHorizAnalysisEnum)); Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 13051) @@ -19,5 +19,5 @@ bool continuous_galerkin=true; int numberofvertices; - bool isstokes,ismacayealpattyn; + bool isstokes,isl1l2,ismacayealpattyn; /*DataSets: */ @@ -27,4 +27,5 @@ iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum); iomodel->Constant(&isstokes,FlowequationIsstokesEnum); + iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); @@ -36,5 +37,5 @@ /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if(!ismacayealpattyn & !isstokes){ + if(!ismacayealpattyn & !isstokes & !isl1l2){ *pnodes=nodes; return; Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/UpdateElementsDiagnosticHoriz.cpp (revision 13051) @@ -19,4 +19,5 @@ int numberofelements; bool ismacayealpattyn; + bool isl1l2; bool isstokes; bool control_analysis; @@ -25,4 +26,5 @@ /*Fetch constants needed: */ iomodel->Constant(&isstokes,FlowequationIsstokesEnum); + iomodel->Constant(&isl1l2,FlowequationIsl1l2Enum); iomodel->Constant(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); iomodel->Constant(&dim,MeshDimensionEnum); @@ -32,5 +34,5 @@ /*Now, is the flag macayaealpattyn on? otherwise, do nothing: */ - if(!ismacayealpattyn & !isstokes) return; + if(!ismacayealpattyn & !isstokes &!isl1l2) return; /*Fetch data needed: */ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 13051) @@ -19,4 +19,7 @@ if (analysis_type==DiagnosticHorizAnalysisEnum){ if (vertices_type[0]==MacAyealApproximationEnum){ + numdofs=2; + } + else if (vertices_type[0]==L1L2ApproximationEnum){ numdofs=2; } Index: /issm/trunk-jpl/src/c/solutions/convergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutions/convergence.cpp (revision 13050) +++ /issm/trunk-jpl/src/c/solutions/convergence.cpp (revision 13051) @@ -81,9 +81,9 @@ //print if(res " << eps_res*100 << " %"); + if(VerboseConvergence()) _pprintLine_(setw(50)< "<parameters->FindParam(&ishutter,FlowequationIshutterEnum); femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); + femmodel->parameters->FindParam(&isl1l2,FlowequationIsl1l2Enum); femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum); femmodel->parameters->FindParam(&isnewton,DiagnosticIsnewtonEnum); @@ -70,5 +72,5 @@ } - if (ismacayealpattyn ^ isstokes){ // ^ = xor + if ((ismacayealpattyn || isl1l2) ^ isstokes){ // ^ = xor if(VerboseSolution()) _pprintLine_(" computing velocities");