Index: /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17113) +++ /issm/trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp (revision 17114) @@ -59,5 +59,22 @@ /*Finite element Analysis*/ void LevelsetAnalysis::Core(FemModel* femmodel){/*{{{*/ - _error_("not implemented yet"); + + /*parameters: */ + bool save_results; + + /*activate formulation: */ + femmodel->SetCurrentConfiguration(LevelsetAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _printf0_("call computational core:\n"); + solutionsequence_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs = MaskIceLevelsetEnum; + femmodel->RequestedOutputsx(&femmodel->results,&outputs,1); + } }/*}}}*/ ElementVector* LevelsetAnalysis::CreateDVector(Element* element){/*{{{*/ @@ -70,10 +87,11 @@ }/*}}}*/ ElementMatrix* LevelsetAnalysis::CreateKMatrix(Element* element){/*{{{*/ + /*Intermediaries */ - - int dim = 2; // solve for LSF in horizontal plane only + int i,j,dim = 2; // solve for LSF in horizontal plane only IssmDouble kappa; IssmDouble Jdet, dt, D_scalar; - IssmDouble u,v,um,vm,ub,vb,vx,vy; + IssmDouble h,hx,hy,hz; + IssmDouble vel,vx,vy,bx,by; IssmDouble* xyz_list = NULL; @@ -86,5 +104,5 @@ IssmDouble* B = xNew(dim*numnodes); IssmDouble* Bprime = xNew(dim*numnodes); - IssmDouble D[dim][dim]; + IssmDouble D[dim][dim], K[dim][dim]; /*Retrieve all inputs and parameters*/ @@ -93,6 +111,4 @@ Input* vx_input = element->GetInput(VxEnum); _assert_(vx_input); Input* vy_input = element->GetInput(VyEnum); _assert_(vy_input); - Input* vxm_input = element->GetInput(VxMeshEnum); _assert_(vxm_input); - Input* vym_input = element->GetInput(VyMeshEnum); _assert_(vym_input); /* Start looping on the number of gaussian points: */ @@ -117,12 +133,12 @@ GetB(B,element,xyz_list,gauss); GetBprime(Bprime,element,xyz_list,gauss); - vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um; - vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm; - ub=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here - vb=0.; - D[0][0]=D_scalar*(vx+ub); + vx_input->GetInputValue(&vx,gauss); // in 3D case, add mesh velocity + vy_input->GetInputValue(&vy,gauss); + bx=0.; // horizontal mass change velocities (melt/refreeze/calving) FIXME: insert values from model here + by=0.; + D[0][0]=D_scalar*(vx+bx); D[0][1]=0.; D[1][0]=0.; - D[1][1]=D_scalar*(vy+vb); + D[1][1]=D_scalar*(vy+by); TripleMultiply(B,dim,numnodes,1, &D[0][0],dim,dim,0, @@ -130,13 +146,45 @@ &Ke->values[0],1); - /* Artificial Diffusion */ - kappa=0.; //FIXME: insert suitable value for kappa - GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed? - D[0][0]=D_scalar*kappa; - D[1][1]=D_scalar*kappa; - TripleMultiply(Bprime,dim,numnodes,1, - &D[0][0],dim,dim,0, - Bprime,dim,numnodes,0, - &Ke->values[0],1); + /* Stabilization */ + int stabilization=1; + switch(stabilization){ + case 0: + // no stabilization, do nothing + break; + case 1: + /* Artificial Diffusion */ + element->ElementSizes(&hx,&hy,&hz); + vel=sqrt(vx*vx + vy*vy ); + h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); //FIXME: is this correct? + + kappa=h*vel/2.; //FIXME: insert suitable value for kappa + //GetBprime(Bprime,element,xyz_list,gauss); // recalculation of Bprime needed? + D[0][0]=D_scalar*kappa; + D[0][1]=0.; + D[1][0]=0.; + D[1][1]=D_scalar*kappa; + TripleMultiply(Bprime,dim,numnodes,1, + &D[0][0],dim,dim,0, + Bprime,dim,numnodes,0, + &Ke->values[0],1); + break; + case 2: + /* Streamline Upwinding */ + element->ElementSizes(&hx,&hy,&hz); + vel=sqrt(vx*vx + vy*vy )+1.e-14; + h=sqrt( pow(hx*vx/vel,2) + pow(hy*vy/vel,2) ); + K[0][0]=h/(2.*vel)*vx*vx; K[0][1]=h/(2.*vel)*vx*vy; + K[1][0]=h/(2.*vel)*vy*vx; K[1][1]=h/(2.*vel)*vy*vy; + for(i=0;ivalues[0],1); + break; + default: + _error_("unknown type of stabilization in LevelsetAnalysis.cpp"); + } } @@ -161,5 +209,5 @@ /*Initialize Element vector*/ - ElementVector* pe = element->NewElementVector(); + ElementVector* pe = element->NewElementVector(); element->FindParam(&dt,TimesteppingTimeStepEnum); @@ -188,5 +236,7 @@ delete gauss; } - else{for(i=0;ivalues[i]=0.;} + else + for(i=0;ivalues[i]=0.; return pe; }/*}}}*/ @@ -195,5 +245,16 @@ }/*}}}*/ void LevelsetAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + + int meshtype; + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DhorizontalEnum: + element->InputUpdateFromSolutionOneDof(solution,MaskIceLevelsetEnum); + break; + case Mesh3DEnum: + element->InputUpdateFromSolutionOneDofCollapsed(solution,MaskIceLevelsetEnum); + break; + default: _error_("mesh "<xmax) xmax=xyz_list[i][0]; + if(xyz_list[i][1]ymax) ymax=xyz_list[i][1]; + } + + *hx=xmax-xmin; + *hy=ymax-ymin; + *hz=0.; } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17113) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17114) @@ -68,5 +68,5 @@ void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); void Delta18oParameterization(void); - void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; + void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz); void ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 17113) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 17114) @@ -50,4 +50,5 @@ femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); @@ -161,7 +162,7 @@ if(islevelset){ if(VerboseSolution()) _printf0_(" computing movement of ice boundaries\n"); - //analysis = new LevelsetAnalysis(); - //analysis->Core(femmodel); - //delete analysis; + analysis = new LevelsetAnalysis(); + analysis->Core(femmodel); + delete analysis; } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 17113) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 17114) @@ -68,4 +68,5 @@ parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum)); parameters->AddObject(iomodel->CopyConstantObject(TransientIsgiaEnum)); + parameters->AddObject(iomodel->CopyConstantObject(TransientIslevelsetEnum)); parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum)); parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum)); Index: /issm/trunk-jpl/src/m/classes/mask.m =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.m (revision 17113) +++ /issm/trunk-jpl/src/m/classes/mask.m (revision 17114) @@ -25,8 +25,8 @@ md = checkfield(md,'fieldname','mask.groundedice_levelset','size',[md.mesh.numberofvertices 1]); md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices 1]); - isice=(md.mask.ice_levelset>0); - if any(sum(isice(md.mesh.elements),2)==0), - warning('elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0'); - end + %isice=(md.mask.ice_levelset>0); + %if any(sum(isice(md.mesh.elements),2)==0), + % warning('elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0'); + % end end % }}} function disp(obj) % {{{ Index: /issm/trunk-jpl/src/m/classes/mask.py =================================================================== --- /issm/trunk-jpl/src/m/classes/mask.py (revision 17113) +++ /issm/trunk-jpl/src/m/classes/mask.py (revision 17114) @@ -33,8 +33,8 @@ md = checkfield(md,'fieldname','mask.ice_levelset' ,'size',[md.mesh.numberofvertices]) - isice=numpy.array(md.mask.ice_levelset>0,int) - totallyicefree=(numpy.sum(isice[md.mesh.elements-1],axis=1)==0).astype(int) - if any(totallyicefree): - raise TypeError("elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0") + #isice=numpy.array(md.mask.ice_levelset<0,int) + #totallyicefree=(numpy.sum(isice[md.mesh.elements-1],axis=1)==0).astype(int) + #if any(totallyicefree): + # raise TypeError("elements with no ice not implemented yet, each element should have at least one vertex with md.mask.ice_levelset > 0") return md Index: /issm/trunk-jpl/src/m/classes/transient.m =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.m (revision 17113) +++ /issm/trunk-jpl/src/m/classes/transient.m (revision 17114) @@ -39,7 +39,9 @@ end % }}} function list = defaultoutputs(self,md) % {{{ - - list = {'SurfaceforcingsMassBalance'}; - + if(self.ismasstransport) + list = {'SurfaceforcingsMassBalance'}; + else + list = {}; + end end % }}} function md = checkconsistency(obj,md,solution,analyses) % {{{ Index: /issm/trunk-jpl/src/m/classes/transient.py =================================================================== --- /issm/trunk-jpl/src/m/classes/transient.py (revision 17113) +++ /issm/trunk-jpl/src/m/classes/transient.py (revision 17114) @@ -38,5 +38,8 @@ def defaultoutputs(self,md): # {{{ - return ['SurfaceforcingsMassBalance'] + if self.ismasstransport: + return ['SurfaceforcingsMassBalance'] + else: + return [] #}}} Index: /issm/trunk-jpl/src/m/parameterization/setmask.m =================================================================== --- /issm/trunk-jpl/src/m/parameterization/setmask.m (revision 17113) +++ /issm/trunk-jpl/src/m/parameterization/setmask.m (revision 17114) @@ -42,5 +42,5 @@ %level sets -md.mask.ice_levelset=ones(md.mesh.numberofvertices,1); +md.mask.ice_levelset=-1.*ones(md.mesh.numberofvertices,1); md.mask.groundedice_levelset=vertexongroundedice; md.mask.groundedice_levelset(find(vertexongroundedice==0.))=-1.; Index: /issm/trunk-jpl/src/m/parameterization/setmask.py =================================================================== --- /issm/trunk-jpl/src/m/parameterization/setmask.py (revision 17113) +++ /issm/trunk-jpl/src/m/parameterization/setmask.py (revision 17114) @@ -42,5 +42,5 @@ #level sets - md.mask.ice_levelset = numpy.ones(md.mesh.numberofvertices,bool) + md.mask.ice_levelset = -1.*numpy.ones((md.mesh.numberofvertices,1)) md.mask.groundedice_levelset = -1.*numpy.ones((md.mesh.numberofvertices,1)) md.mask.groundedice_levelset[md.mesh.elements[numpy.nonzero(elementongroundedice),:]-1]=1.