Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 19385) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp (revision 19386) @@ -115,5 +115,5 @@ /* no source term; damage handled in stress balance */ - f[i]=0; + f[i]=0.; } @@ -632,5 +632,5 @@ }/*}}}*/ void DamageEvolutionAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ - _error_("not implemented yet"); + element->GetSolutionFromInputsOneDof(solution,DamageDbarEnum); }/*}}}*/ void DamageEvolutionAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_index){/*{{{*/ @@ -686,2 +686,172 @@ return; }/*}}}*/ + +/*Flux Correction Transport*/ +ElementMatrix* DamageEvolutionAnalysis::CreateFctKMatrix(Element* element){/*{{{*/ + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Intermediaries */ + IssmDouble Jdet; + IssmDouble vx,vy; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + int dim = 2; + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Ke = element->NewElementMatrix(); + IssmDouble* B = xNew(dim*numnodes); + IssmDouble* Bprime = xNew(dim*numnodes); + IssmDouble* D = xNewZeroInit(dim*dim); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + Input* vxaverage_input=element->GetInput(VxEnum); _assert_(vxaverage_input); + Input* vyaverage_input=element->GetInput(VyEnum); _assert_(vyaverage_input); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + GetB(B,element,dim,xyz_list,gauss); + GetBprime(Bprime,element,dim,xyz_list,gauss); + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + + D[0*dim+0] = -gauss->weight*vx*Jdet; + D[1*dim+1] = -gauss->weight*vy*Jdet; + + TripleMultiply(B,dim,numnodes,1, + D,dim,dim,0, + Bprime,dim,numnodes,0, + &Ke->values[0],1); + + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(B); + xDelete(Bprime); + xDelete(D); + delete gauss; + return Ke; +}/*}}}*/ +ElementMatrix* DamageEvolutionAnalysis::CreateMassMatrix(Element* element){/*{{{*/ + + /* Check if ice in element */ + if(!element->IsIceInElement()) return NULL; + + /*Intermediaries*/ + IssmDouble D,Jdet; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize Element vector and other vectors*/ + ElementMatrix* Me = element->NewElementMatrix(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + + D=gauss->weight*Jdet; + TripleMultiply(basis,1,numnodes,1, + &D,1,1,0, + basis,1,numnodes,0, + &Me->values[0],1); + } + + /*Clean up and return*/ + xDelete(xyz_list); + xDelete(basis); + delete gauss; + return Me; +}/*}}}*/ +void DamageEvolutionAnalysis::FctKMatrix(Matrix** pKff,Matrix** pKfs,FemModel* femmodel){/*{{{*/ + + /*Output*/ + Matrix* Kff = NULL; + Matrix* Kfs = NULL; + + /*Initialize Jacobian Matrix*/ + AllocateSystemMatricesx(&Kff,&Kfs,NULL,NULL,femmodel); + + /*Create and assemble matrix*/ + for(int i=0;ielements->Size();i++){ + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + ElementMatrix* Ke = this->CreateFctKMatrix(element); + if(Ke) Ke->AddToGlobal(Kff,Kfs); + delete Ke; + } + Kff->Assemble(); + Kfs->Assemble(); + + /*Assign output pointer*/ + *pKff=Kff; + if(pKfs){ + *pKfs=Kfs; + } + else{ + delete Kfs; + } +}/*}}}*/ +void DamageEvolutionAnalysis::LumpedMassMatrix(Vector** pMlff,FemModel* femmodel){/*{{{*/ + + /*Intermediaries*/ + int configuration_type; + + /*Initialize Lumped mass matrix (actually we just save its diagonal)*/ + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + int fsize = femmodel->nodes->NumberOfDofs(configuration_type,FsetEnum); + int flocalsize = femmodel->nodes->NumberOfDofsLocal(configuration_type,FsetEnum); + Vector* Mlff = new Vector(flocalsize,fsize); + + /*Create and assemble matrix*/ + for(int i=0;ielements->Size();i++){ + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + ElementMatrix* MLe = this->CreateMassMatrix(element); + if(MLe){ + MLe->Lump(); + MLe->AddDiagonalToGlobal(Mlff); + } + delete MLe; + } + Mlff->Assemble(); + + /*Assign output pointer*/ + *pMlff=Mlff; +}/*}}}*/ +void DamageEvolutionAnalysis::MassMatrix(Matrix** pMff,FemModel* femmodel){/*{{{*/ + + /*Initialize Mass matrix*/ + Matrix *Mff = NULL; + AllocateSystemMatricesx(&Mff,NULL,NULL,NULL,femmodel); + + /*Create and assemble matrix*/ + for(int i=0;ielements->Size();i++){ + Element* element = xDynamicCast(femmodel->elements->GetObjectByOffset(i)); + ElementMatrix* MLe = this->CreateMassMatrix(element); + if(MLe){ + MLe->AddToGlobal(Mff); + } + delete MLe; + } + Mff->Assemble(); + + /*Assign output pointer*/ + *pMff=Mff; +}/*}}}*/ Index: /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 19385) +++ /issm/trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.h (revision 19386) @@ -35,4 +35,11 @@ void InputUpdateFromSolution(IssmDouble* solution,Element* element); void UpdateConstraints(FemModel* femmodel); + + /*FCT*/ + ElementMatrix* CreateFctKMatrix(Element* element); + ElementMatrix* CreateMassMatrix(Element* element); + void FctKMatrix(Matrix** pKff,Matrix** pKfs,FemModel* femmodel); + void LumpedMassMatrix(Vector** pMLff,FemModel* femmodel); + void MassMatrix(Matrix** pMff,FemModel* femmodel); }; #endif Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 19385) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 19386) @@ -876,5 +876,5 @@ iomodel->Constant(&materials_type,MaterialsEnum); if(materials_type==MatdamageiceEnum){ - parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum)); + parameters->AddObject(iomodel->CopyConstantObject(DamageLawEnum)); parameters->AddObject(iomodel->CopyConstantObject(DamageKappaEnum)); parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum)); Index: /issm/trunk-jpl/src/c/cores/damage_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 19385) +++ /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 19386) @@ -1,4 +1,4 @@ /* - * \brief: damgage_core.cpp: core for the damage solution + * \brief: damage_core.cpp: core for the damage solution */ @@ -15,5 +15,5 @@ bool save_results; bool dakota_analysis = false; - int solution_type; + int solution_type,stabilization; int numoutputs = 0; char **requested_outputs = NULL; @@ -25,8 +25,14 @@ femmodel->parameters->FindParam(&numoutputs,DamageEvolutionNumRequestedOutputsEnum); if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DamageEvolutionRequestedOutputsEnum); + femmodel->parameters->FindParam(&stabilization,DamageStabilizationEnum); if(VerboseSolution()) _printf0_(" computing damage\n"); femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum); - solutionsequence_linear(femmodel); + if(stabilization==4){ + solutionsequence_fct(femmodel); + } + else{ + solutionsequence_linear(femmodel); + } if(save_results){ Index: /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp =================================================================== --- /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 19385) +++ /issm/trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp (revision 19386) @@ -1,4 +1,4 @@ -/*!\file: solutionsequence_linear.cpp - * \brief: numerical core of linear solutions +/*!\file: solutionsequence_fct.cpp + * \brief: numerical core of flux corrected transport solution */ @@ -7,4 +7,5 @@ #include "../shared/shared.h" #include "../modules/modules.h" +#include "../analyses/analyses.h" #ifdef _HAVE_PETSC_ @@ -119,5 +120,5 @@ }/*}}}*/ void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/ - /*Create Left Hand side of Lower order solution + /*Create Right Hand side of Lower order solution * * RHS = [ML + (1 − theta) deltaT L^n] u^n @@ -166,5 +167,5 @@ }/*}}}*/ void RichardsonUdot(Vec* pudot,Vec u,Vec Ml,Mat K,Mat Mc){/*{{{*/ - /*Use Richardson's formulato get udot using 5 steps and an initial guess of 0 + /*Use Richardson's formula to get udot using 5 steps and an initial guess of 0 * * udot_new = udot_old + Ml^-1 (K^(n+1) u - Mc udot_old) @@ -336,9 +337,9 @@ }/*}}}*/ #endif -void solutionsequence_fct(FemModel* femmodel){ +void solutionsequence_fct(FemModel* femmodel){/*{{{*/ /*intermediary: */ IssmDouble theta,deltat,dmax; - int configuration_type; + int configuration_type,analysis_type; Vector* Ml = NULL; Matrix* K = NULL; @@ -346,7 +347,6 @@ Vector* ug = NULL; Vector* uf = NULL; - - /*Create analysis*/ - MasstransportAnalysis* analysis = new MasstransportAnalysis(); + MasstransportAnalysis* manalysis = NULL; + DamageEvolutionAnalysis* danalysis = NULL; /*Recover parameters: */ @@ -356,9 +356,23 @@ theta = 0.5; - /*Create lumped mass matrix*/ - analysis->LumpedMassMatrix(&Ml,femmodel); - analysis->MassMatrix(&Mc,femmodel); - analysis->FctKMatrix(&K,NULL,femmodel); - delete analysis; + /*Create analysis*/ + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + switch(analysis_type){ + case MasstransportAnalysisEnum: + manalysis = new MasstransportAnalysis(); + manalysis->LumpedMassMatrix(&Ml,femmodel); + manalysis->MassMatrix(&Mc,femmodel); + manalysis->FctKMatrix(&K,NULL,femmodel); + break; + case DamageEvolutionAnalysisEnum: + danalysis = new DamageEvolutionAnalysis(); + danalysis->LumpedMassMatrix(&Ml,femmodel); + danalysis->MassMatrix(&Mc,femmodel); + danalysis->FctKMatrix(&K,NULL,femmodel); + break; + default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n"); + } + delete manalysis; + delete danalysis; #ifdef _HAVE_PETSC_ @@ -442,3 +456,3 @@ _error_("PETSc needs to be installed"); #endif -} +}/*}}}*/ Index: /issm/trunk-jpl/src/m/classes/damage.m =================================================================== --- /issm/trunk-jpl/src/m/classes/damage.m (revision 19385) +++ /issm/trunk-jpl/src/m/classes/damage.m (revision 19386) @@ -142,5 +142,5 @@ md = checkfield(md,'fieldname','damage.spcdamage','timeseries',1); md = checkfield(md,'fieldname','damage.max_damage','<',1,'>=',0); - md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2]); + md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0 1 2 4]); md = checkfield(md,'fieldname','damage.maxiter','>=0',0); md = checkfield(md,'fieldname','damage.elementinterp','values',{'P1','P2'}); @@ -178,5 +178,5 @@ fielddisplay(self,'max_damage','maximum possible damage (0<=max_damage<1)'); - fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG'); + fielddisplay(self,'stabilization','0: no, 1: artificial_diffusivity, 2: SUPG (not working), 4: flux corrected transport'); fielddisplay(self,'maxiter','maximum number of non linear iterations'); fielddisplay(self,'elementinterp','interpolation scheme for finite elements {''P1'',''P2''}'); Index: /issm/trunk-jpl/src/m/classes/damage.py =================================================================== --- /issm/trunk-jpl/src/m/classes/damage.py (revision 19385) +++ /issm/trunk-jpl/src/m/classes/damage.py (revision 19386) @@ -56,5 +56,5 @@ s+="%s\n" % fielddisplay(self,"max_damage","maximum possible damage (0<=max_damage<1)") - s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG") + s+="%s\n" % fielddisplay(self,"stabilization","0: no, 1: artificial_diffusivity, 2: SUPG (not working), 4: Flux corrected transport") s+="%s\n" % fielddisplay(self,"maxiter","maximum number of non linear iterations") s+="%s\n" % fielddisplay(self,"elementinterp","interpolation scheme for finite elements [''P1'',''P2'']") @@ -126,5 +126,5 @@ md = checkfield(md,'fieldname','damage.law','numel',[1],'values',[0,1,2,3]) md = checkfield(md,'fieldname','damage.spcdamage','timeseries',1) - md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2]) + md = checkfield(md,'fieldname','damage.stabilization','numel',[1],'values',[0,1,2,4]) md = checkfield(md,'fieldname','damage.maxiter','>=0',0) md = checkfield(md,'fieldname','damage.elementinterp','values',['P1','P2']) Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.m =================================================================== --- /issm/trunk-jpl/src/m/mech/mechanicalproperties.m (revision 19385) +++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.m (revision 19386) @@ -62,4 +62,5 @@ nu(location)=B_bar(location); elseif isa(md.materials,'matdamageice') + % FIXME this needs to use md.damage.D or damage from a solution Zinv=md.materials.rheology_Z(index)*summation/3; location=find(second_inv~=0); Index: /issm/trunk-jpl/src/m/mech/mechanicalproperties.py =================================================================== --- /issm/trunk-jpl/src/m/mech/mechanicalproperties.py (revision 19385) +++ /issm/trunk-jpl/src/m/mech/mechanicalproperties.py (revision 19386) @@ -77,4 +77,5 @@ elif 'matdamageice' in md.materials.__module__: print 'computing damage-dependent properties!' + # FIXME might be using damage from a solution, not the initial md.damage.D Zinv=npy.dot(1-md.damage.D[index-1],summation/3.).reshape(-1,) location=npy.nonzero(second_inv)