Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.cpp (revision 18926) @@ -40,43 +40,20 @@ } /*}}}*/ -void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - - switch(this->law){ - case 1: - GetAlpha2Viscous(palpha2,gauss); - break; - case 2: - GetAlpha2Weertman(palpha2,gauss); - break; - case 3: - GetAlpha2Hydro(palpha2,gauss); - break; - case 4: - GetAlpha2Temp(palpha2,gauss); - break; - case 5: - GetAlpha2WaterLayer(palpha2,gauss); - break; - case 6: - GetAlpha2WeertmanTemp(palpha2,gauss); - break; - default: - _error_("not supported"); - } - -}/*}}}*/ -void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - - /*This routine calculates the basal friction coefficient - alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ +void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ + + /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. + * FrictionGetAlphaComplement is used in control methods on drag, and it computes: + * alpha_complement= Neff ^r * vel ^s*/ + + if(this->law!=1)_error_("not supported"); /*diverse: */ IssmDouble r,s; - IssmDouble drag_p, drag_q; + IssmDouble vx,vy,vz,vmag; + IssmDouble drag_p,drag_q; IssmDouble Neff; - IssmDouble thickness,bed; - IssmDouble vx,vy,vz,vmag; IssmDouble drag_coefficient; - IssmDouble alpha2; + IssmDouble bed,thickness; + IssmDouble alpha_complement; /*Recover parameters: */ @@ -98,4 +75,5 @@ if(Neff<0)Neff=0; + //We need the velocity magnitude to evaluate the basal stress: switch(dim){ case 1: @@ -119,51 +97,36 @@ /*Check to prevent dividing by zero if vmag==0*/ - if(vmag==0. && (s-1.)<0.) alpha2=0.; - else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.)); - _assert_(!xIsNan(alpha2)); - - /*Assign output pointers:*/ - *palpha2=alpha2; -}/*}}}*/ -void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ - - /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */ - - /*diverse: */ - IssmDouble C,m; - IssmDouble vx,vy,vz,vmag; - IssmDouble alpha2; - - /*Recover parameters: */ - element->GetInputValue(&C,gauss,FrictionCEnum); - element->GetInputValue(&m,FrictionMEnum); - - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } - - /*Check to prevent dividing by zero if vmag==0*/ - if(vmag==0. && (1./m-1.)<0.) alpha2=0.; - else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.)); - _assert_(!xIsNan(alpha2)); - - /*Assign output pointers:*/ - *palpha2=alpha2; + if(vmag==0. && (s-1.)<0.) alpha_complement=0.; + else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan(alpha_complement)); + + /*Assign output pointers:*/ + *palpha_complement=alpha_complement; +} +/*}}}*/ +void Friction::GetAlpha2(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + switch(this->law){ + case 1: + GetAlpha2Viscous(palpha2,gauss); + break; + case 2: + GetAlpha2Weertman(palpha2,gauss); + break; + case 3: + GetAlpha2Hydro(palpha2,gauss); + break; + case 4: + GetAlpha2Temp(palpha2,gauss); + break; + case 5: + GetAlpha2WaterLayer(palpha2,gauss); + break; + case 6: + GetAlpha2WeertmanTemp(palpha2,gauss); + break; + default: + _error_("not supported"); + } + }/*}}}*/ void Friction::GetAlpha2Hydro(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ @@ -256,4 +219,64 @@ element->parameters->FindParam(&gamma,FrictionGammaEnum); alpha2 = alpha2 / exp((T-Tpmp)/gamma); + + /*Assign output pointers:*/ + *palpha2=alpha2; +}/*}}}*/ +void Friction::GetAlpha2Viscous(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + /*This routine calculates the basal friction coefficient + alpha2= drag^2 * Neff ^r * | vel | ^(s-1), with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p**/ + + /*diverse: */ + IssmDouble r,s; + IssmDouble drag_p, drag_q; + IssmDouble Neff; + IssmDouble thickness,bed; + IssmDouble vx,vy,vz,vmag; + IssmDouble drag_coefficient; + IssmDouble alpha2; + + /*Recover parameters: */ + element->GetInputValue(&drag_p,FrictionPEnum); + element->GetInputValue(&drag_q,FrictionQEnum); + element->GetInputValue(&thickness, gauss,ThicknessEnum); + element->GetInputValue(&bed, gauss,BaseEnum); + element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + + //compute r and q coefficients: */ + r=drag_q/drag_p; + s=1./drag_p; + + //From bed and thickness, compute effective pressure when drag is viscous: + Neff=gravity*(rho_ice*thickness+rho_water*bed); + if(Neff<0)Neff=0; + + switch(dim){ + case 1: + element->GetInputValue(&vx,gauss,VxEnum); + vmag=sqrt(vx*vx); + break; + case 2: + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + vmag=sqrt(vx*vx+vy*vy); + break; + case 3: + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + element->GetInputValue(&vz,gauss,VzEnum); + vmag=sqrt(vx*vx+vy*vy+vz*vz); + break; + default: + _error_("not supported"); + } + + /*Check to prevent dividing by zero if vmag==0*/ + if(vmag==0. && (s-1.)<0.) alpha2=0.; + else alpha2=drag_coefficient*drag_coefficient*pow(Neff,r)*pow(vmag,(s-1.)); + _assert_(!xIsNan(alpha2)); /*Assign output pointers:*/ @@ -322,4 +345,45 @@ *palpha2=alpha2; }/*}}}*/ +void Friction::GetAlpha2Weertman(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ + + /*This routine calculates the basal friction coefficient alpha2= C^-1/m |v|^(1/m-1) */ + + /*diverse: */ + IssmDouble C,m; + IssmDouble vx,vy,vz,vmag; + IssmDouble alpha2; + + /*Recover parameters: */ + element->GetInputValue(&C,gauss,FrictionCEnum); + element->GetInputValue(&m,FrictionMEnum); + + switch(dim){ + case 1: + element->GetInputValue(&vx,gauss,VxEnum); + vmag=sqrt(vx*vx); + break; + case 2: + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + vmag=sqrt(vx*vx+vy*vy); + break; + case 3: + element->GetInputValue(&vx,gauss,VxEnum); + element->GetInputValue(&vy,gauss,VyEnum); + element->GetInputValue(&vz,gauss,VzEnum); + vmag=sqrt(vx*vx+vy*vy+vz*vz); + break; + default: + _error_("not supported"); + } + + /*Check to prevent dividing by zero if vmag==0*/ + if(vmag==0. && (1./m-1.)<0.) alpha2=0.; + else alpha2=pow(C,-1./m)*pow(vmag,(1./m-1.)); + _assert_(!xIsNan(alpha2)); + + /*Assign output pointers:*/ + *palpha2=alpha2; +}/*}}}*/ void Friction::GetAlpha2WeertmanTemp(IssmDouble* palpha2, Gauss* gauss){/*{{{*/ /*Here, we want to parameterize the friction as a function of temperature @@ -349,66 +413,2 @@ *palpha2=alpha2; }/*}}}*/ -void Friction::GetAlphaComplement(IssmDouble* palpha_complement, Gauss* gauss){/*{{{*/ - - /* FrictionGetAlpha2 computes alpha2= drag^2 * Neff ^r * vel ^s, with Neff=rho_ice*g*thickness+rho_ice*g*bed, r=q/p and s=1/p. - * FrictionGetAlphaComplement is used in control methods on drag, and it computes: - * alpha_complement= Neff ^r * vel ^s*/ - - if(this->law!=1)_error_("not supported"); - - /*diverse: */ - IssmDouble r,s; - IssmDouble vx,vy,vz,vmag; - IssmDouble drag_p,drag_q; - IssmDouble Neff; - IssmDouble drag_coefficient; - IssmDouble bed,thickness; - IssmDouble alpha_complement; - - /*Recover parameters: */ - element->GetInputValue(&drag_p,FrictionPEnum); - element->GetInputValue(&drag_q,FrictionQEnum); - element->GetInputValue(&thickness, gauss,ThicknessEnum); - element->GetInputValue(&bed, gauss,BaseEnum); - element->GetInputValue(&drag_coefficient, gauss,FrictionCoefficientEnum); - IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum); - IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); - IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); - - //compute r and q coefficients: */ - r=drag_q/drag_p; - s=1./drag_p; - - //From bed and thickness, compute effective pressure when drag is viscous: - Neff=gravity*(rho_ice*thickness+rho_water*bed); - if(Neff<0)Neff=0; - - //We need the velocity magnitude to evaluate the basal stress: - switch(dim){ - case 1: - element->GetInputValue(&vx,gauss,VxEnum); - vmag=sqrt(vx*vx); - break; - case 2: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - vmag=sqrt(vx*vx+vy*vy); - break; - case 3: - element->GetInputValue(&vx,gauss,VxEnum); - element->GetInputValue(&vy,gauss,VyEnum); - element->GetInputValue(&vz,gauss,VzEnum); - vmag=sqrt(vx*vx+vy*vy+vz*vz); - break; - default: - _error_("not supported"); - } - - /*Check to prevent dividing by zero if vmag==0*/ - if(vmag==0. && (s-1.)<0.) alpha_complement=0.; - else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));_assert_(!xIsNan(alpha_complement)); - - /*Assign output pointers:*/ - *palpha_complement=alpha_complement; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Friction.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Friction.h (revision 18926) @@ -29,12 +29,12 @@ void Echo(void); + void GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss); void GetAlpha2(IssmDouble* palpha2,Gauss* gauss); - void GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss); - void GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Hydro(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2Temp(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2Viscous(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2WaterLayer(IssmDouble* palpha2,Gauss* gauss); + void GetAlpha2Weertman(IssmDouble* palpha2,Gauss* gauss); void GetAlpha2WeertmanTemp(IssmDouble* palpha2,Gauss* gauss); - void GetAlphaComplement(IssmDouble* alpha_complement,Gauss* gauss); }; Index: /issm/trunk-jpl/src/c/classes/Loads/Load.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Load.h (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Load.h (revision 18926) @@ -26,17 +26,17 @@ virtual ~Load(){}; virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; - virtual void ResetHooks()=0; - virtual bool IsPenalty(void)=0; - virtual int GetNumberOfNodes(void)=0; - virtual void GetNodesSidList(int* sidlist)=0; - virtual void GetNodesLidList(int* lidlist)=0; - virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; + virtual void CreateJacobianMatrix(Matrix* Jff)=0; virtual void CreateKMatrix(Matrix* Kff, Matrix* Kfs)=0; virtual void CreatePVector(Vector* pf)=0; - virtual void CreateJacobianMatrix(Matrix* Jff)=0; + virtual void GetNodesLidList(int* lidlist)=0; + virtual void GetNodesSidList(int* sidlist)=0; + virtual int GetNumberOfNodes(void)=0; + virtual bool InAnalysis(int analysis_type)=0; + virtual bool IsPenalty(void)=0; virtual void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax)=0; virtual void PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs, IssmDouble kmax)=0; virtual void PenaltyCreatePVector(Vector* pf, IssmDouble kmax)=0; - virtual bool InAnalysis(int analysis_type)=0; + virtual void ResetHooks()=0; + virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; virtual void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum)=0; }; Index: /issm/trunk-jpl/src/c/classes/Loads/Loads.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Loads.cpp (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Loads.cpp (revision 18926) @@ -48,18 +48,4 @@ } /*}}}*/ -void Loads::ResetHooks(){/*{{{*/ - - vector::iterator object; - Load* load=NULL; - - for ( object=objects.begin() ; object < objects.end(); object++ ){ - - load=xDynamicCast((*object)); - load->ResetHooks(); - - } - -} -/*}}}*/ bool Loads::IsPenalty(int analysis_type){/*{{{*/ @@ -86,5 +72,5 @@ } /*}}}*/ -int Loads::MaxNumNodes(int analysis_type){/*{{{*/ +int Loads::MaxNumNodes(int analysis_type){/*{{{*/ int max=0; @@ -109,5 +95,5 @@ } /*}}}*/ -int Loads::NumberOfLoads(void){/*{{{*/ +int Loads::NumberOfLoads(void){/*{{{*/ int localloads; @@ -124,5 +110,5 @@ } /*}}}*/ -int Loads::NumberOfLoads(int analysis_type){/*{{{*/ +int Loads::NumberOfLoads(int analysis_type){/*{{{*/ int localloads = 0; @@ -145,23 +131,16 @@ } /*}}}*/ -int Loads::Size(void){/*{{{*/ +void Loads::ResetHooks(){/*{{{*/ - return this->DataSet::Size(); -} -/*}}}*/ -int Loads::Size(int analysis_type){/*{{{*/ + vector::iterator object; + Load* load=NULL; - int localloads = 0; + for ( object=objects.begin() ; object < objects.end(); object++ ){ - /*Get number of local loads*/ - for(int i=0;iSize();i++){ + load=xDynamicCast((*object)); + load->ResetHooks(); - Load* load=xDynamicCast(this->GetObjectByOffset(i)); - - /*Check that this load corresponds to our analysis currently being carried out: */ - if (load->InAnalysis(analysis_type)) localloads++; } - return localloads; } /*}}}*/ @@ -180,2 +159,23 @@ } /*}}}*/ +int Loads::Size(void){/*{{{*/ + + return this->DataSet::Size(); +} +/*}}}*/ +int Loads::Size(int analysis_type){/*{{{*/ + + int localloads = 0; + + /*Get number of local loads*/ + for(int i=0;iSize();i++){ + + Load* load=xDynamicCast(this->GetObjectByOffset(i)); + + /*Check that this load corresponds to our analysis currently being carried out: */ + if (load->InAnalysis(analysis_type)) localloads++; + } + + return localloads; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Loads.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Loads.h (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Loads.h (revision 18926) @@ -24,9 +24,9 @@ /*numerics*/ void Configure(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); - void ResetHooks(); bool IsPenalty(int analysis); int MaxNumNodes(int analysis); int NumberOfLoads(void); int NumberOfLoads(int analysis); + void ResetHooks(); void SetCurrentConfiguration(Elements* elements,Loads* loads, Nodes* nodes, Vertices* vertices, Materials* materials,Parameters* parameters); int Size(int analysis); Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.cpp (revision 18926) @@ -131,16 +131,32 @@ /*Object virtual functions definitions:*/ -void Numericalflux::Echo(void){/*{{{*/ - _printf_("Numericalflux:\n"); - _printf_(" id: " << id << "\n"); - _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); - _printf_(" flux_type: " << this->flux_type<< "\n"); - hnodes->Echo(); - hvertices->Echo(); - helement->Echo(); - _printf_(" parameters: " << parameters << "\n"); -} -/*}}}*/ -void Numericalflux::DeepEcho(void){/*{{{*/ +Object* Numericalflux::copy() {/*{{{*/ + + Numericalflux* numericalflux=NULL; + + numericalflux=new Numericalflux(); + + /*copy fields: */ + numericalflux->id=this->id; + numericalflux->analysis_type=this->analysis_type; + numericalflux->flux_type=this->flux_type; + + /*point parameters: */ + numericalflux->parameters=this->parameters; + + /*now deal with hooks and objects: */ + numericalflux->hnodes = (Hook*)this->hnodes->copy(); + numericalflux->hvertices = (Hook*)this->hvertices->copy(); + numericalflux->helement = (Hook*)this->helement->copy(); + + /*corresponding fields*/ + numericalflux->nodes = (Node**)numericalflux->hnodes->deliverp(); + numericalflux->vertices = (Vertex**)numericalflux->hvertices->deliverp(); + numericalflux->element = (Element*)numericalflux->helement->delivers(); + + return numericalflux; +} +/*}}}*/ +void Numericalflux::DeepEcho(void){/*{{{*/ _printf_("Numericalflux:\n"); @@ -158,39 +174,23 @@ } /*}}}*/ -int Numericalflux::Id(void){/*{{{*/ +void Numericalflux::Echo(void){/*{{{*/ + _printf_("Numericalflux:\n"); + _printf_(" id: " << id << "\n"); + _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); + _printf_(" flux_type: " << this->flux_type<< "\n"); + hnodes->Echo(); + hvertices->Echo(); + helement->Echo(); + _printf_(" parameters: " << parameters << "\n"); +} +/*}}}*/ +int Numericalflux::Id(void){/*{{{*/ return id; } /*}}}*/ -int Numericalflux::ObjectEnum(void){/*{{{*/ +int Numericalflux::ObjectEnum(void){/*{{{*/ return NumericalfluxEnum; -} -/*}}}*/ -Object* Numericalflux::copy() {/*{{{*/ - - Numericalflux* numericalflux=NULL; - - numericalflux=new Numericalflux(); - - /*copy fields: */ - numericalflux->id=this->id; - numericalflux->analysis_type=this->analysis_type; - numericalflux->flux_type=this->flux_type; - - /*point parameters: */ - numericalflux->parameters=this->parameters; - - /*now deal with hooks and objects: */ - numericalflux->hnodes = (Hook*)this->hnodes->copy(); - numericalflux->hvertices = (Hook*)this->hvertices->copy(); - numericalflux->helement = (Hook*)this->helement->copy(); - - /*corresponding fields*/ - numericalflux->nodes = (Node**)numericalflux->hnodes->deliverp(); - numericalflux->vertices = (Vertex**)numericalflux->hvertices->deliverp(); - numericalflux->element = (Element*)numericalflux->helement->delivers(); - - return numericalflux; } /*}}}*/ @@ -212,22 +212,4 @@ /*point parameters to real dataset: */ this->parameters=parametersin; -} -/*}}}*/ -void Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ - -} -/*}}}*/ -void Numericalflux::ResetHooks(){/*{{{*/ - - this->nodes=NULL; - this->vertices=NULL; - this->element=NULL; - this->parameters=NULL; - - /*Get Element type*/ - this->hnodes->reset(); - this->hvertices->reset(); - this->helement->reset(); - } /*}}}*/ @@ -291,4 +273,21 @@ } /*}}}*/ +void Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/ + + _assert_(lidlist); + _assert_(nodes); + + switch(this->flux_type){ + case InternalEnum: + for(int i=0;iLid(); + return; + case BoundaryEnum: + for(int i=0;iLid(); + return; + default: + _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); + } +} +/*}}}*/ void Numericalflux::GetNodesSidList(int* sidlist){/*{{{*/ @@ -308,21 +307,4 @@ } /*}}}*/ -void Numericalflux::GetNodesLidList(int* lidlist){/*{{{*/ - - _assert_(lidlist); - _assert_(nodes); - - switch(this->flux_type){ - case InternalEnum: - for(int i=0;iLid(); - return; - case BoundaryEnum: - for(int i=0;iLid(); - return; - default: - _error_("Numericalflux type " << EnumToStringx(this->flux_type) << " not supported yet"); - } -} -/*}}}*/ int Numericalflux::GetNumberOfNodes(void){/*{{{*/ @@ -338,25 +320,43 @@ } /*}}}*/ -bool Numericalflux::IsPenalty(void){/*{{{*/ - return false; -} -/*}}}*/ -void Numericalflux::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs,IssmDouble kmax){/*{{{*/ - - /*No stiffness loads applied, do nothing: */ - return; - -} -/*}}}*/ -void Numericalflux::PenaltyCreatePVector(Vector* pf,IssmDouble kmax){/*{{{*/ - - /*No penalty loads applied, do nothing: */ - return; - -} -/*}}}*/ bool Numericalflux::InAnalysis(int in_analysis_type){/*{{{*/ if (in_analysis_type==this->analysis_type) return true; else return false; +} +/*}}}*/ +bool Numericalflux::IsPenalty(void){/*{{{*/ + return false; +} +/*}}}*/ +void Numericalflux::PenaltyCreateKMatrix(Matrix* Kff, Matrix* Kfs,IssmDouble kmax){/*{{{*/ + + /*No stiffness loads applied, do nothing: */ + return; + +} +/*}}}*/ +void Numericalflux::PenaltyCreatePVector(Vector* pf,IssmDouble kmax){/*{{{*/ + + /*No penalty loads applied, do nothing: */ + return; + +} +/*}}}*/ +void Numericalflux::ResetHooks(){/*{{{*/ + + this->nodes=NULL; + this->vertices=NULL; + this->element=NULL; + this->parameters=NULL; + + /*Get Element type*/ + this->hnodes->reset(); + this->hvertices->reset(); + this->helement->reset(); + +} +/*}}}*/ +void Numericalflux::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ + } /*}}}*/ @@ -417,4 +417,175 @@ /*Numericalflux management*/ +ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/ + + switch(this->flux_type){ + case InternalEnum: + return CreateKMatrixAdjointBalancethicknessInternal(); + case BoundaryEnum: + return CreateKMatrixAdjointBalancethicknessBoundary(); + default: + _error_("type not supported yet"); + } +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/ + + ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary(); + if(Ke) Ke->Transpose(); + return Ke; +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/ + + ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal(); + if (Ke) Ke->Transpose(); + return Ke; +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/ + + switch(this->flux_type){ + case InternalEnum: + return CreateKMatrixBalancethicknessInternal(); + case BoundaryEnum: + return CreateKMatrixBalancethicknessBoundary(); + default: + _error_("type not supported yet"); + } +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/ + + /* constants*/ + const int numdof=NDOF1*NUMNODES_BOUNDARY; + + /* Intermediaries*/ + int i,j,ig,index1,index2; + IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble L[numdof]; + IssmDouble Ke_g[numdof][numdof]; + GaussTria *gauss; + + /*Initialize Element matrix and return if necessary*/ + ElementMatrix* Ke = NULL; + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* vxaverage_input=tria->inputs->GetInput(VxEnum); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); + GetNormal(&normal[0],xyz_list); + + /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + + gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + vxaverage_input->GetInputValue(&mean_vx,gauss); + vyaverage_input->GetInputValue(&mean_vy,gauss); + delete gauss; + + UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if (UdotN<=0){ + return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ + } + else{ + Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters); + } + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + UdotN=vx*normal[0]+vy*normal[1]; + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL=gauss->weight*Jdet*UdotN; + + TripleMultiply(&L[0],1,numdof,1, + &DL,1,1,0, + &L[0],1,numdof,0, + &Ke_g[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/ + + /* constants*/ + const int numdof=NDOF1*NUMNODES_INTERNAL; + + /* Intermediaries*/ + int i,j,ig,index1,index2; + IssmDouble DL1,DL2,Jdet,vx,vy,UdotN; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble B[numdof]; + IssmDouble Bprime[numdof]; + IssmDouble Ke_g1[numdof][numdof]; + IssmDouble Ke_g2[numdof][numdof]; + GaussTria *gauss; + + /*Initialize Element matrix and return if necessary*/ + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + Input* vxaverage_input=tria->inputs->GetInput(VxEnum); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); + GetNormal(&normal[0],xyz_list); + + /* Start looping on the number of gaussian points: */ + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement()); + tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement()); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + UdotN=vx*normal[0]+vy*normal[1]; + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL1=gauss->weight*Jdet*UdotN/2; + DL2=gauss->weight*Jdet*fabs(UdotN)/2; + + TripleMultiply(&B[0],1,numdof,1, + &DL1,1,1,0, + &Bprime[0],1,numdof,0, + &Ke_g1[0][0],0); + TripleMultiply(&B[0],1,numdof,1, + &DL2,1,1,0, + &B[0],1,numdof,0, + &Ke_g2[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_g1[i][j]; + for(i=0;ivalues[i*numdof+j]+=Ke_g2[i][j]; + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ ElementMatrix* Numericalflux::CreateKMatrixMasstransport(void){/*{{{*/ @@ -427,4 +598,75 @@ _error_("type not supported yet"); } +} +/*}}}*/ +ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/ + + /* constants*/ + const int numdof=NDOF1*NUMNODES_BOUNDARY; + + /* Intermediaries*/ + int i,j,ig,index1,index2; + IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble L[numdof]; + IssmDouble Ke_g[numdof][numdof]; + GaussTria *gauss; + + /*Initialize Element matrix and return if necessary*/ + ElementMatrix* Ke = NULL; + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); + Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); + GetNormal(&normal[0],xyz_list); + + /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + + gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + vxaverage_input->GetInputValue(&mean_vx,gauss); + vyaverage_input->GetInputValue(&mean_vy,gauss); + delete gauss; + + UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if (UdotN<=0){ + return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ + } + else{ + Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters); + } + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + UdotN=vx*normal[0]+vy*normal[1]; + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL=gauss->weight*Jdet*dt*UdotN; + + TripleMultiply(&L[0],1,numdof,1, + &DL,1,1,0, + &L[0],1,numdof,0, + &Ke_g[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; + } + + /*Clean up and return*/ + delete gauss; + return Ke; } /*}}}*/ @@ -493,332 +735,8 @@ } /*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixMasstransportBoundary(void){/*{{{*/ - - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; - - /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble L[numdof]; - IssmDouble Ke_g[numdof][numdof]; - GaussTria *gauss; - - /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke = NULL; - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vxaverage_input=tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); - Input* vyaverage_input=tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); - GetNormal(&normal[0],xyz_list); - - /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); - gauss->GaussEdgeCenter(index1,index2); - vxaverage_input->GetInputValue(&mean_vx,gauss); - vyaverage_input->GetInputValue(&mean_vy,gauss); - delete gauss; - - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN<=0){ - return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ - } - else{ - Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - UdotN=vx*normal[0]+vy*normal[1]; - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - DL=gauss->weight*Jdet*dt*UdotN; - - TripleMultiply(&L[0],1,numdof,1, - &DL,1,1,0, - &L[0],1,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixBalancethickness(void){/*{{{*/ - - switch(this->flux_type){ - case InternalEnum: - return CreateKMatrixBalancethicknessInternal(); - case BoundaryEnum: - return CreateKMatrixBalancethicknessBoundary(); - default: - _error_("type not supported yet"); - } -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessInternal(void){/*{{{*/ - - /* constants*/ - const int numdof=NDOF1*NUMNODES_INTERNAL; - - /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL1,DL2,Jdet,vx,vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble B[numdof]; - IssmDouble Bprime[numdof]; - IssmDouble Ke_g1[numdof][numdof]; - IssmDouble Ke_g2[numdof][numdof]; - GaussTria *gauss; - - /*Initialize Element matrix and return if necessary*/ - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - ElementMatrix* Ke=new ElementMatrix(nodes,NUMNODES_INTERNAL,this->parameters); - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vxaverage_input=tria->inputs->GetInput(VxEnum); - Input* vyaverage_input=tria->inputs->GetInput(VyEnum); - GetNormal(&normal[0],xyz_list); - - /* Start looping on the number of gaussian points: */ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetSegmentBFlux(&B[0],gauss,index1,index2,tria->FiniteElement()); - tria->GetSegmentBprimeFlux(&Bprime[0],gauss,index1,index2,tria->FiniteElement()); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - UdotN=vx*normal[0]+vy*normal[1]; - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - DL1=gauss->weight*Jdet*UdotN/2; - DL2=gauss->weight*Jdet*fabs(UdotN)/2; - - TripleMultiply(&B[0],1,numdof,1, - &DL1,1,1,0, - &Bprime[0],1,numdof,0, - &Ke_g1[0][0],0); - TripleMultiply(&B[0],1,numdof,1, - &DL2,1,1,0, - &B[0],1,numdof,0, - &Ke_g2[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g1[i][j]; - for(i=0;ivalues[i*numdof+j]+=Ke_g2[i][j]; - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixBalancethicknessBoundary(void){/*{{{*/ - - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; - - /* Intermediaries*/ - int i,j,ig,index1,index2; - IssmDouble DL,Jdet,vx,vy,mean_vx,mean_vy,UdotN; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble L[numdof]; - IssmDouble Ke_g[numdof][numdof]; - GaussTria *gauss; - - /*Initialize Element matrix and return if necessary*/ - ElementMatrix* Ke = NULL; - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - Input* vxaverage_input=tria->inputs->GetInput(VxEnum); - Input* vyaverage_input=tria->inputs->GetInput(VyEnum); - GetNormal(&normal[0],xyz_list); - - /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); - gauss->GaussEdgeCenter(index1,index2); - vxaverage_input->GetInputValue(&mean_vx,gauss); - vyaverage_input->GetInputValue(&mean_vy,gauss); - delete gauss; - - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN<=0){ - return NULL; /*(u,n)<0 -> inflow, PenaltyCreatePVector will take care of it*/ - } - else{ - Ke=new ElementMatrix(nodes,NUMNODES_BOUNDARY,this->parameters); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - UdotN=vx*normal[0]+vy*normal[1]; - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - DL=gauss->weight*Jdet*UdotN; - - TripleMultiply(&L[0],1,numdof,1, - &DL,1,1,0, - &L[0],1,numdof,0, - &Ke_g[0][0],0); - - for(i=0;ivalues[i*numdof+j]+=Ke_g[i][j]; - } - - /*Clean up and return*/ - delete gauss; - return Ke; -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethickness(void){/*{{{*/ - - switch(this->flux_type){ - case InternalEnum: - return CreateKMatrixAdjointBalancethicknessInternal(); - case BoundaryEnum: - return CreateKMatrixAdjointBalancethicknessBoundary(); - default: - _error_("type not supported yet"); - } -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessInternal(void){/*{{{*/ - - ElementMatrix* Ke=CreateKMatrixBalancethicknessInternal(); - if (Ke) Ke->Transpose(); - return Ke; -} -/*}}}*/ -ElementMatrix* Numericalflux::CreateKMatrixAdjointBalancethicknessBoundary(void){/*{{{*/ - - ElementMatrix* Ke=CreateKMatrixBalancethicknessBoundary(); - if(Ke) Ke->Transpose(); - return Ke; -} -/*}}}*/ -ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/ - - switch(this->flux_type){ - case InternalEnum: - return CreatePVectorMasstransportInternal(); - case BoundaryEnum: - return CreatePVectorMasstransportBoundary(); - default: - _error_("type not supported yet"); - } -} -/*}}}*/ -ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/ - - /*Nothing added to PVector*/ +ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/ + + /*No PVector for the Adjoint*/ return NULL; - -} -/*}}}*/ -ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/ - - /* constants*/ - const int numdof=NDOF1*NUMNODES_BOUNDARY; - - /* Intermediaries*/ - int i,ig,index1,index2; - IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness; - IssmDouble xyz_list[NUMVERTICES][3]; - IssmDouble normal[2]; - IssmDouble L[numdof]; - GaussTria *gauss; - - /*Initialize Load Vector and return if necessary*/ - ElementVector* pe = NULL; - Tria* tria=(Tria*)element; - if(!tria->IsIceInElement()) return NULL; - - /*Retrieve all inputs and parameters*/ - GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); - parameters->FindParam(&dt,TimesteppingTimeStepEnum); - Input* vxaverage_input =tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); - Input* vyaverage_input =tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); - Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input); - GetNormal(&normal[0],xyz_list); - - /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ - index1=tria->GetNodeIndex(nodes[0]); - index2=tria->GetNodeIndex(nodes[1]); - - gauss=new GaussTria(); - gauss->GaussEdgeCenter(index1,index2); - vxaverage_input->GetInputValue(&mean_vx,gauss); - vyaverage_input->GetInputValue(&mean_vy,gauss); - delete gauss; - - UdotN=mean_vx*normal[0]+mean_vy*normal[1]; - if (UdotN>0){ - return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ - } - else{ - pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters); - } - - /* Start looping on the number of gaussian points: */ - gauss=new GaussTria(index1,index2,2); - for(ig=gauss->begin();igend();ig++){ - - gauss->GaussPoint(ig); - - tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); - - vxaverage_input->GetInputValue(&vx,gauss); - vyaverage_input->GetInputValue(&vy,gauss); - spcthickness_input->GetInputValue(&thickness,gauss); - if(xIsNan(thickness)) _error_("Cannot weakly apply constraint because NaN was provided"); - - UdotN=vx*normal[0]+vy*normal[1]; - tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); - DL= - gauss->weight*Jdet*dt*UdotN*thickness; - - for(i=0;ivalues[i] += DL*L[i]; - } - - /*Clean up and return*/ - delete gauss; - return pe; } /*}}}*/ @@ -833,11 +751,4 @@ _error_("type not supported yet"); } -} -/*}}}*/ -ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/ - - /*Nothing added to PVector*/ - return NULL; - } /*}}}*/ @@ -908,11 +819,100 @@ } /*}}}*/ -ElementVector* Numericalflux::CreatePVectorAdjointBalancethickness(void){/*{{{*/ - - /*No PVector for the Adjoint*/ +ElementVector* Numericalflux::CreatePVectorBalancethicknessInternal(void){/*{{{*/ + + /*Nothing added to PVector*/ return NULL; -} -/*}}}*/ -void Numericalflux:: GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/ + +} +/*}}}*/ +ElementVector* Numericalflux::CreatePVectorMasstransport(void){/*{{{*/ + + switch(this->flux_type){ + case InternalEnum: + return CreatePVectorMasstransportInternal(); + case BoundaryEnum: + return CreatePVectorMasstransportBoundary(); + default: + _error_("type not supported yet"); + } +} +/*}}}*/ +ElementVector* Numericalflux::CreatePVectorMasstransportBoundary(void){/*{{{*/ + + /* constants*/ + const int numdof=NDOF1*NUMNODES_BOUNDARY; + + /* Intermediaries*/ + int i,ig,index1,index2; + IssmDouble DL,Jdet,dt,vx,vy,mean_vx,mean_vy,UdotN,thickness; + IssmDouble xyz_list[NUMVERTICES][3]; + IssmDouble normal[2]; + IssmDouble L[numdof]; + GaussTria *gauss; + + /*Initialize Load Vector and return if necessary*/ + ElementVector* pe = NULL; + Tria* tria=(Tria*)element; + if(!tria->IsIceInElement()) return NULL; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + parameters->FindParam(&dt,TimesteppingTimeStepEnum); + Input* vxaverage_input =tria->inputs->GetInput(VxEnum); _assert_(vxaverage_input); + Input* vyaverage_input =tria->inputs->GetInput(VyEnum); _assert_(vyaverage_input); + Input* spcthickness_input=tria->inputs->GetInput(MasstransportSpcthicknessEnum); _assert_(spcthickness_input); + GetNormal(&normal[0],xyz_list); + + /*Check wether it is an inflow or outflow BC (0 is the middle of the segment)*/ + index1=tria->GetNodeIndex(nodes[0]); + index2=tria->GetNodeIndex(nodes[1]); + + gauss=new GaussTria(); + gauss->GaussEdgeCenter(index1,index2); + vxaverage_input->GetInputValue(&mean_vx,gauss); + vyaverage_input->GetInputValue(&mean_vy,gauss); + delete gauss; + + UdotN=mean_vx*normal[0]+mean_vy*normal[1]; + if (UdotN>0){ + return NULL; /*(u,n)>0 -> outflow, PenaltyCreateKMatrix will take care of it*/ + } + else{ + pe=new ElementVector(nodes,NUMNODES_BOUNDARY,this->parameters); + } + + /* Start looping on the number of gaussian points: */ + gauss=new GaussTria(index1,index2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + tria->GetSegmentNodalFunctions(&L[0],gauss,index1,index2,tria->FiniteElement()); + + vxaverage_input->GetInputValue(&vx,gauss); + vyaverage_input->GetInputValue(&vy,gauss); + spcthickness_input->GetInputValue(&thickness,gauss); + if(xIsNan(thickness)) _error_("Cannot weakly apply constraint because NaN was provided"); + + UdotN=vx*normal[0]+vy*normal[1]; + tria->GetSegmentJacobianDeterminant(&Jdet,&xyz_list[0][0],gauss); + DL= - gauss->weight*Jdet*dt*UdotN*thickness; + + for(i=0;ivalues[i] += DL*L[i]; + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +ElementVector* Numericalflux::CreatePVectorMasstransportInternal(void){/*{{{*/ + + /*Nothing added to PVector*/ + return NULL; + +} +/*}}}*/ +void Numericalflux::GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]){/*{{{*/ /*Build unit outward pointing vector*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Numericalflux.h (revision 18926) @@ -40,54 +40,54 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ + Object *copy(); + void DeepEcho(); void Echo(); - void DeepEcho(); int Id(); int ObjectEnum(); - Object *copy(); /*}}}*/ /*Update virtual functions resolution: {{{*/ - void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/} - void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/} - void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/} void InputUpdateFromConstant(IssmDouble constant, int name){/*Do nothing*/}; void InputUpdateFromConstant(int constant, int name){/*Do nothing*/}; void InputUpdateFromConstant(bool constant, int name){_error_("Not implemented yet!");} void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*Do nothing*/} + void InputUpdateFromVector(IssmDouble* vector, int name, int type){/*Do nothing*/} + void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*Do nothing*/} /*}}}*/ /*Load virtual functions definitions: {{{*/ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void ResetHooks(); + void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; void CreateKMatrix(Matrix* Kff, Matrix* Kfs); void CreatePVector(Vector* pf); + void GetNodesLidList(int* lidlist); void GetNodesSidList(int* sidlist); - void GetNodesLidList(int* lidlist); int GetNumberOfNodes(void); - void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; + bool InAnalysis(int analysis_type); bool IsPenalty(void); void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){_error_("Not implemented yet");}; void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); + void ResetHooks(); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); - bool InAnalysis(int analysis_type); + void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); /*}}}*/ /*Numericalflux management:{{{*/ + ElementMatrix* CreateKMatrixAdjointBalancethickness(void); + ElementMatrix* CreateKMatrixAdjointBalancethicknessBoundary(void); + ElementMatrix* CreateKMatrixAdjointBalancethicknessInternal(void); + ElementMatrix* CreateKMatrixBalancethickness(void); + ElementMatrix* CreateKMatrixBalancethicknessBoundary(void); + ElementMatrix* CreateKMatrixBalancethicknessInternal(void); + ElementMatrix* CreateKMatrixMasstransport(void); + ElementMatrix* CreateKMatrixMasstransportBoundary(void); + ElementMatrix* CreateKMatrixMasstransportInternal(void); + ElementVector* CreatePVectorAdjointBalancethickness(void); + ElementVector* CreatePVectorBalancethickness(void); + ElementVector* CreatePVectorBalancethicknessBoundary(void); + ElementVector* CreatePVectorBalancethicknessInternal(void); + ElementVector* CreatePVectorMasstransport(void); + ElementVector* CreatePVectorMasstransportBoundary(void); + ElementVector* CreatePVectorMasstransportInternal(void); void GetNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]); - ElementMatrix* CreateKMatrixMasstransport(void); - ElementMatrix* CreateKMatrixMasstransportInternal(void); - ElementMatrix* CreateKMatrixMasstransportBoundary(void); - ElementMatrix* CreateKMatrixBalancethickness(void); - ElementMatrix* CreateKMatrixBalancethicknessInternal(void); - ElementMatrix* CreateKMatrixBalancethicknessBoundary(void); - ElementMatrix* CreateKMatrixAdjointBalancethickness(void); - ElementMatrix* CreateKMatrixAdjointBalancethicknessInternal(void); - ElementMatrix* CreateKMatrixAdjointBalancethicknessBoundary(void); - ElementVector* CreatePVectorMasstransport(void); - ElementVector* CreatePVectorMasstransportInternal(void); - ElementVector* CreatePVectorMasstransportBoundary(void); - ElementVector* CreatePVectorBalancethickness(void); - ElementVector* CreatePVectorBalancethicknessInternal(void); - ElementVector* CreatePVectorBalancethicknessBoundary(void); - ElementVector* CreatePVectorAdjointBalancethickness(void); /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Pengrid.cpp (revision 18926) @@ -81,9 +81,36 @@ /*Object virtual functions definitions:*/ -void Pengrid::Echo(void){/*{{{*/ - this->DeepEcho(); -} -/*}}}*/ -void Pengrid::DeepEcho(void){/*{{{*/ +Object* Pengrid::copy() {/*{{{*/ + + Pengrid* pengrid=NULL; + + pengrid=new Pengrid(); + + /*copy fields: */ + pengrid->id=this->id; + pengrid->analysis_type=this->analysis_type; + + /*point parameters: */ + pengrid->parameters=this->parameters; + + /*now deal with hooks and objects: */ + pengrid->hnode=(Hook*)this->hnode->copy(); + pengrid->hmatpar=(Hook*)this->hmatpar->copy(); + pengrid->helement=(Hook*)this->helement->copy(); + + /*corresponding fields*/ + pengrid->node =(Node*)pengrid->hnode->delivers(); + pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers(); + pengrid->element=(Element*)pengrid->helement->delivers(); + + //let's not forget internals + pengrid->active=this->active=0; + pengrid->zigzag_counter=this->zigzag_counter=0; + + return pengrid; + +} +/*}}}*/ +void Pengrid::DeepEcho(void){/*{{{*/ _printf_("Pengrid:\n"); @@ -99,40 +126,13 @@ } /*}}}*/ -int Pengrid::Id(void){ return id; }/*{{{*/ -/*}}}*/ -int Pengrid::ObjectEnum(void){/*{{{*/ +void Pengrid::Echo(void){/*{{{*/ + this->DeepEcho(); +} +/*}}}*/ +int Pengrid::Id(void){ return id; }/*{{{*/ +/*}}}*/ +int Pengrid::ObjectEnum(void){/*{{{*/ return PengridEnum; -} -/*}}}*/ -Object* Pengrid::copy() {/*{{{*/ - - Pengrid* pengrid=NULL; - - pengrid=new Pengrid(); - - /*copy fields: */ - pengrid->id=this->id; - pengrid->analysis_type=this->analysis_type; - - /*point parameters: */ - pengrid->parameters=this->parameters; - - /*now deal with hooks and objects: */ - pengrid->hnode=(Hook*)this->hnode->copy(); - pengrid->hmatpar=(Hook*)this->hmatpar->copy(); - pengrid->helement=(Hook*)this->helement->copy(); - - /*corresponding fields*/ - pengrid->node =(Node*)pengrid->hnode->delivers(); - pengrid->matpar =(Matpar*)pengrid->hmatpar->delivers(); - pengrid->element=(Element*)pengrid->helement->delivers(); - - //let's not forget internals - pengrid->active=this->active=0; - pengrid->zigzag_counter=this->zigzag_counter=0; - - return pengrid; - } /*}}}*/ @@ -154,22 +154,4 @@ /*point parameters to real dataset: */ this->parameters=parametersin; -} -/*}}}*/ -void Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ - -} -/*}}}*/ -void Pengrid::ResetHooks(){/*{{{*/ - - this->node=NULL; - this->element=NULL; - this->matpar=NULL; - this->parameters=NULL; - - /*Get Element type*/ - this->hnode->reset(); - this->helement->reset(); - this->hmatpar->reset(); - } /*}}}*/ @@ -203,4 +185,12 @@ } /*}}}*/ +void Pengrid::GetNodesLidList(int* lidlist){/*{{{*/ + + _assert_(lidlist); + _assert_(node); + + lidlist[0]=node->Lid(); +} +/*}}}*/ void Pengrid::GetNodesSidList(int* sidlist){/*{{{*/ @@ -211,15 +201,16 @@ } /*}}}*/ -void Pengrid::GetNodesLidList(int* lidlist){/*{{{*/ - - _assert_(lidlist); - _assert_(node); - - lidlist[0]=node->Lid(); -} -/*}}}*/ int Pengrid::GetNumberOfNodes(void){/*{{{*/ return NUMVERTICES; +} +/*}}}*/ +bool Pengrid::InAnalysis(int in_analysis_type){/*{{{*/ + if (in_analysis_type==this->analysis_type)return true; + else return false; +} +/*}}}*/ +bool Pengrid::IsPenalty(void){/*{{{*/ + return true; } /*}}}*/ @@ -282,11 +273,20 @@ } /*}}}*/ -bool Pengrid::InAnalysis(int in_analysis_type){/*{{{*/ - if (in_analysis_type==this->analysis_type)return true; - else return false; -} -/*}}}*/ -bool Pengrid::IsPenalty(void){/*{{{*/ - return true; +void Pengrid::ResetHooks(){/*{{{*/ + + this->node=NULL; + this->element=NULL; + this->matpar=NULL; + this->parameters=NULL; + + /*Get Element type*/ + this->hnode->reset(); + this->helement->reset(); + this->hmatpar->reset(); + +} +/*}}}*/ +void Pengrid::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ + } /*}}}*/ @@ -343,16 +343,4 @@ /*Update virtual functions definitions:*/ -void Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/ - /*Nothing updated yet*/ -} -/*}}}*/ -void Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ - /*Nothing updated yet*/ -} -/*}}}*/ -void Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/ - /*Nothing updated yet*/ -} -/*}}}*/ void Pengrid::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/ /*Nothing*/ @@ -374,7 +362,19 @@ } /*}}}*/ +void Pengrid::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/ + /*Nothing updated yet*/ +} +/*}}}*/ +void Pengrid::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/ + /*Nothing updated yet*/ +} +/*}}}*/ +void Pengrid::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/ + /*Nothing updated yet*/ +} +/*}}}*/ /*Pengrid management:*/ -void Pengrid::ConstraintActivate(int* punstable){/*{{{*/ +void Pengrid::ConstraintActivate(int* punstable){/*{{{*/ int analysis_type; @@ -404,179 +404,5 @@ } /*}}}*/ -void Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/ - - // The penalty is stable if it doesn't change during to successive iterations. - IssmDouble pressure; - IssmDouble temperature; - IssmDouble t_pmp; - int new_active; - int unstable=0; - int penalty_lock; - - /*recover pointers: */ - Penta* penta=(Penta*)element; - - /*check that pengrid is not a clone (penalty to be added only once)*/ - if (node->IsClone()){ - unstable=0; - *punstable=unstable; - return; - } - - //First recover pressure and temperature values, using the element: */ - penta->GetInputValue(&pressure,node,PressureEnum); - penta->GetInputValue(&temperature,node,TemperaturePicardEnum); - - //Recover our data: - parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum); - - //Compute pressure melting point - t_pmp=matpar->TMeltingPoint(pressure); - - //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated. - - if (temperature>t_pmp){ - new_active=1; - } - else{ - new_active=0; - } - - //Figure out stability of this penalty - if (active==new_active){ - unstable=0; - } - else{ - unstable=1; - if(penalty_lock)zigzag_counter++; - } - - /*If penalty keeps zigzagging more than 5 times: */ - if(penalty_lock){ - if(zigzag_counter>penalty_lock){ - unstable=0; - active=1; - } - } - - //Set penalty flag - active=new_active; - - //*Assign output pointers:*/ - *punstable=unstable; -} -/*}}}*/ -ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){/*{{{*/ - - IssmDouble pressure,temperature,t_pmp; - IssmDouble penalty_factor; - - Penta* penta=(Penta*)element; - - /*check that pengrid is not a clone (penalty to be added only once)*/ - if (node->IsClone()) return NULL; - ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters); - - /*Retrieve all parameters*/ - penta->GetInputValue(&pressure,node,PressureEnum); - penta->GetInputValue(&temperature,node,TemperatureEnum); - parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); - - /*Compute pressure melting point*/ - t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; - - /*Add penalty load*/ - if (temperaturevalues[0]=kmax*pow(10.,penalty_factor); - } - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ -ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(IssmDouble kmax){/*{{{*/ - - IssmDouble penalty_factor; - - /*Initialize Element matrix and return if necessary*/ - if(!this->active) return NULL; - ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters); - - /*recover parameters: */ - parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); - - Ke->values[0]=kmax*pow(10.,penalty_factor); - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ -ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){/*{{{*/ - - IssmDouble pressure; - IssmDouble temperature; - IssmDouble melting_offset; - IssmDouble t_pmp; - IssmDouble dt,penalty_factor; - - /*recover pointers: */ - Penta* penta=(Penta*)element; - - /*check that pengrid is not a clone (penalty to be added only once)*/ - if (node->IsClone()) return NULL; - ElementVector* pe=new ElementVector(&node,NUMVERTICES,this->parameters); - - /*Retrieve all parameters*/ - penta->GetInputValue(&pressure,node,PressureEnum); - penta->GetInputValue(&temperature,node,TemperatureEnum); - parameters->FindParam(&melting_offset,MeltingOffsetEnum); - parameters->FindParam(&dt,TimesteppingTimeStepEnum); - parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); - - /*Compute pressure melting point*/ - t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; - - /*Add penalty load - This time, the penalty must have the same value as the one used for the thermal computation - so that the corresponding melting can be computed correctly - In the thermal computation, we used kmax=melting_offset, and the same penalty_factor*/ - if (temperaturevalues[0]=0; - } - else{ - if (reCast(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt; - else pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp); - } - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -ElementVector* Pengrid::PenaltyCreatePVectorThermal(IssmDouble kmax){/*{{{*/ - - IssmDouble pressure; - IssmDouble t_pmp; - IssmDouble penalty_factor; - - Penta* penta=(Penta*)element; - - /*Initialize Element matrix and return if necessary*/ - if(!this->active) return NULL; - ElementVector* pe=new ElementVector(&node,1,this->parameters); - - /*Retrieve all parameters*/ - penta->GetInputValue(&pressure,node,PressureEnum); - parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); - - /*Compute pressure melting point*/ - t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; - - pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp; - - /*Clean up and return*/ - return pe; -} -/*}}}*/ -void Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/ +void Pengrid::ConstraintActivateHydrologyDCInefficient(int* punstable){/*{{{*/ // The penalty is stable if it doesn't change during two consecutive iterations. @@ -638,4 +464,83 @@ } /*}}}*/ +void Pengrid::ConstraintActivateThermal(int* punstable){/*{{{*/ + + // The penalty is stable if it doesn't change during to successive iterations. + IssmDouble pressure; + IssmDouble temperature; + IssmDouble t_pmp; + int new_active; + int unstable=0; + int penalty_lock; + + /*recover pointers: */ + Penta* penta=(Penta*)element; + + /*check that pengrid is not a clone (penalty to be added only once)*/ + if (node->IsClone()){ + unstable=0; + *punstable=unstable; + return; + } + + //First recover pressure and temperature values, using the element: */ + penta->GetInputValue(&pressure,node,PressureEnum); + penta->GetInputValue(&temperature,node,TemperaturePicardEnum); + + //Recover our data: + parameters->FindParam(&penalty_lock,ThermalPenaltyLockEnum); + + //Compute pressure melting point + t_pmp=matpar->TMeltingPoint(pressure); + + //Figure out if temperature is over melting_point, in which case, this penalty needs to be activated. + + if (temperature>t_pmp){ + new_active=1; + } + else{ + new_active=0; + } + + //Figure out stability of this penalty + if (active==new_active){ + unstable=0; + } + else{ + unstable=1; + if(penalty_lock)zigzag_counter++; + } + + /*If penalty keeps zigzagging more than 5 times: */ + if(penalty_lock){ + if(zigzag_counter>penalty_lock){ + unstable=0; + active=1; + } + } + + //Set penalty flag + active=new_active; + + //*Assign output pointers:*/ + *punstable=unstable; +} +/*}}}*/ +ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/ + + IssmDouble moulin_load,dt; + + /*Initialize Element matrix*/ + ElementVector* pe=new ElementVector(&node,1,this->parameters); + + this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum); + parameters->FindParam(&dt,TimesteppingTimeStepEnum); + + if(dt!=0.0) pe->values[0]=moulin_load*dt; + + /*Clean up and return*/ + return pe; + } +/*}}}*/ ElementMatrix* Pengrid::PenaltyCreateKMatrixHydrologyDCInefficient(IssmDouble kmax){/*{{{*/ IssmDouble penalty_factor; @@ -647,4 +552,49 @@ if(!this->active) return NULL; ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters); + + Ke->values[0]=kmax*pow(10.,penalty_factor); + + /*Clean up and return*/ + return Ke; +} +/*}}}*/ +ElementMatrix* Pengrid::PenaltyCreateKMatrixMelting(IssmDouble kmax){/*{{{*/ + + IssmDouble pressure,temperature,t_pmp; + IssmDouble penalty_factor; + + Penta* penta=(Penta*)element; + + /*check that pengrid is not a clone (penalty to be added only once)*/ + if (node->IsClone()) return NULL; + ElementMatrix* Ke=new ElementMatrix(&node,1,this->parameters); + + /*Retrieve all parameters*/ + penta->GetInputValue(&pressure,node,PressureEnum); + penta->GetInputValue(&temperature,node,TemperatureEnum); + parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); + + /*Compute pressure melting point*/ + t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + + /*Add penalty load*/ + if (temperaturevalues[0]=kmax*pow(10.,penalty_factor); + } + + /*Clean up and return*/ + return Ke; +} +/*}}}*/ +ElementMatrix* Pengrid::PenaltyCreateKMatrixThermal(IssmDouble kmax){/*{{{*/ + + IssmDouble penalty_factor; + + /*Initialize Element matrix and return if necessary*/ + if(!this->active) return NULL; + ElementMatrix* Ke=new ElementMatrix(&node,NUMVERTICES,this->parameters); + + /*recover parameters: */ + parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); Ke->values[0]=kmax*pow(10.,penalty_factor); @@ -678,21 +628,71 @@ } /*}}}*/ -ElementVector* Pengrid::CreatePVectorHydrologyDCInefficient(void){/*{{{*/ - - IssmDouble moulin_load,dt; - - /*Initialize Element matrix*/ - ElementVector* pe=new ElementVector(&node,1,this->parameters); - - this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum); +ElementVector* Pengrid::PenaltyCreatePVectorMelting(IssmDouble kmax){/*{{{*/ + + IssmDouble pressure; + IssmDouble temperature; + IssmDouble melting_offset; + IssmDouble t_pmp; + IssmDouble dt,penalty_factor; + + /*recover pointers: */ + Penta* penta=(Penta*)element; + + /*check that pengrid is not a clone (penalty to be added only once)*/ + if (node->IsClone()) return NULL; + ElementVector* pe=new ElementVector(&node,NUMVERTICES,this->parameters); + + /*Retrieve all parameters*/ + penta->GetInputValue(&pressure,node,PressureEnum); + penta->GetInputValue(&temperature,node,TemperatureEnum); + parameters->FindParam(&melting_offset,MeltingOffsetEnum); parameters->FindParam(&dt,TimesteppingTimeStepEnum); - - if(dt!=0.0) pe->values[0]=moulin_load*dt; + parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); + + /*Compute pressure melting point*/ + t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + + /*Add penalty load + This time, the penalty must have the same value as the one used for the thermal computation + so that the corresponding melting can be computed correctly + In the thermal computation, we used kmax=melting_offset, and the same penalty_factor*/ + if (temperaturevalues[0]=0; + } + else{ + if (reCast(dt)) pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp)/dt; + else pe->values[0]=melting_offset*pow(10.,penalty_factor)*(temperature-t_pmp); + } /*Clean up and return*/ return pe; - } -/*}}}*/ -void Pengrid::ResetConstraint(void){/*{{{*/ +} +/*}}}*/ +ElementVector* Pengrid::PenaltyCreatePVectorThermal(IssmDouble kmax){/*{{{*/ + + IssmDouble pressure; + IssmDouble t_pmp; + IssmDouble penalty_factor; + + Penta* penta=(Penta*)element; + + /*Initialize Element matrix and return if necessary*/ + if(!this->active) return NULL; + ElementVector* pe=new ElementVector(&node,1,this->parameters); + + /*Retrieve all parameters*/ + penta->GetInputValue(&pressure,node,PressureEnum); + parameters->FindParam(&penalty_factor,ThermalPenaltyFactorEnum); + + /*Compute pressure melting point*/ + t_pmp=matpar->GetMeltingPoint()-matpar->GetBeta()*pressure; + + pe->values[0]=kmax*pow(10.,penalty_factor)*t_pmp; + + /*Clean up and return*/ + return pe; +} +/*}}}*/ +void Pengrid::ResetConstraint(void){/*{{{*/ active = 0; zigzag_counter = 0; Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.cpp (revision 18926) @@ -45,5 +45,36 @@ /*Object virtual functions definitions:*/ -void Penpair::Echo(void){/*{{{*/ +Object* Penpair::copy() {/*{{{*/ + + Penpair* penpair=NULL; + + penpair=new Penpair(); + + /*copy fields: */ + penpair->id=this->id; + penpair->analysis_type=this->analysis_type; + + /*now deal with hooks and objects: */ + penpair->hnodes=(Hook*)this->hnodes->copy(); + penpair->nodes =(Node**)penpair->hnodes->deliverp(); + + /*point parameters: */ + penpair->parameters=this->parameters; + + return penpair; + +} +/*}}}*/ +void Penpair::DeepEcho(void){/*{{{*/ + + _printf_("Penpair:\n"); + _printf_(" id: " << id << "\n"); + _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); + hnodes->DeepEcho(); + + return; +} +/*}}}*/ +void Penpair::Echo(void){/*{{{*/ _printf_("Penpair:\n"); @@ -55,40 +86,9 @@ } /*}}}*/ -void Penpair::DeepEcho(void){/*{{{*/ - - _printf_("Penpair:\n"); - _printf_(" id: " << id << "\n"); - _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); - hnodes->DeepEcho(); - - return; -} -/*}}}*/ -int Penpair::Id(void){ return id; }/*{{{*/ -/*}}}*/ -int Penpair::ObjectEnum(void){/*{{{*/ +int Penpair::Id(void){ return id; }/*{{{*/ +/*}}}*/ +int Penpair::ObjectEnum(void){/*{{{*/ return PenpairEnum; -} -/*}}}*/ -Object* Penpair::copy() {/*{{{*/ - - Penpair* penpair=NULL; - - penpair=new Penpair(); - - /*copy fields: */ - penpair->id=this->id; - penpair->analysis_type=this->analysis_type; - - /*now deal with hooks and objects: */ - penpair->hnodes=(Hook*)this->hnodes->copy(); - penpair->nodes =(Node**)penpair->hnodes->deliverp(); - - /*point parameters: */ - penpair->parameters=this->parameters; - - return penpair; - } /*}}}*/ @@ -109,16 +109,6 @@ } /*}}}*/ -void Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ - -} -/*}}}*/ -void Penpair::ResetHooks(){/*{{{*/ - - this->nodes=NULL; - this->parameters=NULL; - - /*Get Element type*/ - this->hnodes->reset(); - +void Penpair::CreateJacobianMatrix(Matrix* Jff){/*{{{*/ + this->CreateKMatrix(Jff,NULL); } /*}}}*/ @@ -137,8 +127,4 @@ } /*}}}*/ -void Penpair::CreateJacobianMatrix(Matrix* Jff){/*{{{*/ - this->CreateKMatrix(Jff,NULL); -} -/*}}}*/ void Penpair::GetNodesSidList(int* sidlist){/*{{{*/ @@ -162,6 +148,15 @@ } /*}}}*/ +bool Penpair::InAnalysis(int in_analysis_type){/*{{{*/ + if (in_analysis_type==this->analysis_type)return true; + else return false; +} +/*}}}*/ bool Penpair::IsPenalty(void){/*{{{*/ return true; +} +/*}}}*/ +void Penpair::PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){/*{{{*/ + this->PenaltyCreateKMatrix(Jff,NULL,kmax); } /*}}}*/ @@ -196,11 +191,16 @@ } /*}}}*/ -void Penpair::PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){/*{{{*/ - this->PenaltyCreateKMatrix(Jff,NULL,kmax); -} -/*}}}*/ -bool Penpair::InAnalysis(int in_analysis_type){/*{{{*/ - if (in_analysis_type==this->analysis_type)return true; - else return false; +void Penpair::ResetHooks(){/*{{{*/ + + this->nodes=NULL; + this->parameters=NULL; + + /*Get Element type*/ + this->hnodes->reset(); + +} +/*}}}*/ +void Penpair::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ + } /*}}}*/ @@ -279,4 +279,56 @@ /*Penpair management:*/ +ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/ + + const int numdof=NUMVERTICES*NDOF1; + IssmDouble penalty_factor; + + /*Initialize Element vector and return if necessary*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters); + + /*recover parameters: */ + parameters->FindParam(&penalty_factor,MasstransportPenaltyFactorEnum); + + //Create elementary matrix: add penalty to + Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor); + Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor); + Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor); + Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor); + + /*Clean up and return*/ + return Ke; +} +/*}}}*/ +ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax){/*{{{*/ + + const int numdof=NUMVERTICES*NDOF3; + IssmDouble penalty_offset; + + /*Initialize Element vector and return if necessary*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters); + + /*recover parameters: */ + parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum); + + //Create elementary matrix: add penalty to + Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset); + Ke->values[0*numdof+3]=-kmax*pow(10.,penalty_offset); + Ke->values[3*numdof+0]=-kmax*pow(10.,penalty_offset); + Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset); + + Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset); + Ke->values[1*numdof+4]=-kmax*pow(10.,penalty_offset); + Ke->values[4*numdof+1]=-kmax*pow(10.,penalty_offset); + Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset); + + Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset); + Ke->values[2*numdof+5]=-kmax*pow(10.,penalty_offset); + Ke->values[5*numdof+2]=-kmax*pow(10.,penalty_offset); + Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset); + + /*Clean up and return*/ + return Ke; +} +/*}}}*/ ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax){/*{{{*/ @@ -338,54 +390,2 @@ } /*}}}*/ -ElementMatrix* Penpair::PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax){/*{{{*/ - - const int numdof=NUMVERTICES*NDOF3; - IssmDouble penalty_offset; - - /*Initialize Element vector and return if necessary*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters); - - /*recover parameters: */ - parameters->FindParam(&penalty_offset,StressbalancePenaltyFactorEnum); - - //Create elementary matrix: add penalty to - Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_offset); - Ke->values[0*numdof+3]=-kmax*pow(10.,penalty_offset); - Ke->values[3*numdof+0]=-kmax*pow(10.,penalty_offset); - Ke->values[3*numdof+3]=+kmax*pow(10.,penalty_offset); - - Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_offset); - Ke->values[1*numdof+4]=-kmax*pow(10.,penalty_offset); - Ke->values[4*numdof+1]=-kmax*pow(10.,penalty_offset); - Ke->values[4*numdof+4]=+kmax*pow(10.,penalty_offset); - - Ke->values[2*numdof+2]=+kmax*pow(10.,penalty_offset); - Ke->values[2*numdof+5]=-kmax*pow(10.,penalty_offset); - Ke->values[5*numdof+2]=-kmax*pow(10.,penalty_offset); - Ke->values[5*numdof+5]=+kmax*pow(10.,penalty_offset); - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ -ElementMatrix* Penpair::PenaltyCreateKMatrixMasstransport(IssmDouble kmax){/*{{{*/ - - const int numdof=NUMVERTICES*NDOF1; - IssmDouble penalty_factor; - - /*Initialize Element vector and return if necessary*/ - ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters); - - /*recover parameters: */ - parameters->FindParam(&penalty_factor,MasstransportPenaltyFactorEnum); - - //Create elementary matrix: add penalty to - Ke->values[0*numdof+0]=+kmax*pow(10.,penalty_factor); - Ke->values[0*numdof+1]=-kmax*pow(10.,penalty_factor); - Ke->values[1*numdof+0]=-kmax*pow(10.,penalty_factor); - Ke->values[1*numdof+1]=+kmax*pow(10.,penalty_factor); - - /*Clean up and return*/ - return Ke; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Penpair.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Penpair.h (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Penpair.h (revision 18926) @@ -31,41 +31,41 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); - int ObjectEnum(); Object* copy(); + void DeepEcho(); + void Echo(); + int Id(); + int ObjectEnum(); /*}}}*/ /*Update virtual functions resolution: {{{*/ - void InputUpdateFromVector(IssmDouble* vector, int name, int type); - void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols,int name, int type){_error_("Not implemented yet!");} - void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");} void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name); void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrow, int ncols,int name, int type){_error_("Not implemented yet!");} + void InputUpdateFromVector(IssmDouble* vector, int name, int type); + void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");} /*}}}*/ /*Load virtual functions definitions: {{{*/ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); + void CreateJacobianMatrix(Matrix* Jff); + void CreateKMatrix(Matrix* Kff, Matrix* Kfs); + void CreatePVector(Vector* pf); + void GetNodesLidList(int* lidlist); + void GetNodesSidList(int* sidlist); + int GetNumberOfNodes(void); + bool InAnalysis(int analysis_type); + bool IsPenalty(void); + void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax); + void PenaltyCreateKMatrix(Matrix* Kff,Matrix* Kfs,IssmDouble kmax); + void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); void ResetHooks(); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void CreateKMatrix(Matrix* Kff, Matrix* Kfs); - void CreatePVector(Vector* pf); - void CreateJacobianMatrix(Matrix* Jff); - void GetNodesSidList(int* sidlist); - void GetNodesLidList(int* lidlist); - int GetNumberOfNodes(void); - bool IsPenalty(void); - void PenaltyCreateKMatrix(Matrix* Kff,Matrix* Kfs,IssmDouble kmax); - void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); - void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); - bool InAnalysis(int analysis_type); /*}}}*/ /*Penpair management: {{{*/ + ElementMatrix* PenaltyCreateKMatrixMasstransport(IssmDouble kmax); + ElementMatrix* PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax); ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax); ElementMatrix* PenaltyCreateKMatrixStressbalanceSSAHO(IssmDouble kmax); - ElementMatrix* PenaltyCreateKMatrixStressbalanceFS(IssmDouble kmax); - ElementMatrix* PenaltyCreateKMatrixMasstransport(IssmDouble kmax); /*}}}*/ }; Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.cpp (revision 18926) @@ -108,5 +108,62 @@ /*Object virtual functions definitions:*/ -void Riftfront::Echo(void){/*{{{*/ +Object* Riftfront::copy() {/*{{{*/ + + Riftfront* riftfront=NULL; + + riftfront=new Riftfront(); + + /*copy fields: */ + riftfront->id=this->id; + riftfront->analysis_type=this->analysis_type; + riftfront->type=this->type; + riftfront->fill=this->fill; + riftfront->friction=this->friction; + riftfront->fractionincrement=this->fractionincrement; + riftfront->shelf=this->shelf; + + /*point parameters: */ + riftfront->parameters=this->parameters; + + /*now deal with hooks and objects: */ + riftfront->hnodes=(Hook*)this->hnodes->copy(); + riftfront->helements=(Hook*)this->helements->copy(); + riftfront->hmatpar=(Hook*)this->hmatpar->copy(); + + /*corresponding fields*/ + riftfront->nodes =(Node**)riftfront->hnodes->deliverp(); + riftfront->elements=(Element**)riftfront->helements->deliverp(); + riftfront->matpar =(Matpar*)riftfront->hmatpar->delivers(); + + /*internal data: */ + riftfront->penalty_lock=this->penalty_lock; + riftfront->active=this->active; + riftfront->frozen=this->frozen; + riftfront->state=this->state; + riftfront->counter=this->counter; + riftfront->prestable=this->prestable; + riftfront->material_converged=this->material_converged; + riftfront->normal[0]=this->normal[0]; + riftfront->normal[1]=this->normal[1]; + riftfront->length=this->length; + riftfront->fraction=this->fraction; + + return riftfront; + +} +/*}}}*/ +void Riftfront::DeepEcho(void){/*{{{*/ + + _printf_("Riftfront:\n"); + _printf_(" id: " << id << "\n"); + _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); + hnodes->DeepEcho(); + helements->DeepEcho(); + hmatpar->DeepEcho(); + _printf_(" parameters\n"); + if(parameters)parameters->DeepEcho(); +} +/*}}}*/ +void Riftfront::Echo(void){/*{{{*/ _printf_("Riftfront:\n"); @@ -134,66 +191,9 @@ } /*}}}*/ -void Riftfront::DeepEcho(void){/*{{{*/ - - _printf_("Riftfront:\n"); - _printf_(" id: " << id << "\n"); - _printf_(" analysis_type: " << EnumToStringx(analysis_type) << "\n"); - hnodes->DeepEcho(); - helements->DeepEcho(); - hmatpar->DeepEcho(); - _printf_(" parameters\n"); - if(parameters)parameters->DeepEcho(); -} -/*}}}*/ -int Riftfront::Id(void){ return id; }/*{{{*/ -/*}}}*/ -int Riftfront::ObjectEnum(void){/*{{{*/ +int Riftfront::Id(void){ return id; }/*{{{*/ +/*}}}*/ +int Riftfront::ObjectEnum(void){/*{{{*/ return RiftfrontEnum; - -} -/*}}}*/ -Object* Riftfront::copy() {/*{{{*/ - - Riftfront* riftfront=NULL; - - riftfront=new Riftfront(); - - /*copy fields: */ - riftfront->id=this->id; - riftfront->analysis_type=this->analysis_type; - riftfront->type=this->type; - riftfront->fill=this->fill; - riftfront->friction=this->friction; - riftfront->fractionincrement=this->fractionincrement; - riftfront->shelf=this->shelf; - - /*point parameters: */ - riftfront->parameters=this->parameters; - - /*now deal with hooks and objects: */ - riftfront->hnodes=(Hook*)this->hnodes->copy(); - riftfront->helements=(Hook*)this->helements->copy(); - riftfront->hmatpar=(Hook*)this->hmatpar->copy(); - - /*corresponding fields*/ - riftfront->nodes =(Node**)riftfront->hnodes->deliverp(); - riftfront->elements=(Element**)riftfront->helements->deliverp(); - riftfront->matpar =(Matpar*)riftfront->hmatpar->delivers(); - - /*internal data: */ - riftfront->penalty_lock=this->penalty_lock; - riftfront->active=this->active; - riftfront->frozen=this->frozen; - riftfront->state=this->state; - riftfront->counter=this->counter; - riftfront->prestable=this->prestable; - riftfront->material_converged=this->material_converged; - riftfront->normal[0]=this->normal[0]; - riftfront->normal[1]=this->normal[1]; - riftfront->length=this->length; - riftfront->fraction=this->fraction; - - return riftfront; } @@ -234,24 +234,42 @@ } /*}}}*/ -void Riftfront::ResetHooks(){/*{{{*/ - - this->nodes=NULL; - this->elements=NULL; - this->matpar=NULL; - this->parameters=NULL; - - /*Get Element type*/ - this->hnodes->reset(); - this->helements->reset(); - this->hmatpar->reset(); - +void Riftfront::CreateKMatrix(Matrix* Kff, Matrix* Kfs){/*{{{*/ + /*do nothing: */ + return; +} +/*}}}*/ +void Riftfront::CreatePVector(Vector* pf){/*{{{*/ + /*do nothing: */ + return; +} +/*}}}*/ +void Riftfront::GetNodesLidList(int* lidlist){/*{{{*/ + + _assert_(lidlist); + _assert_(nodes); + + for(int i=0;iLid(); +} +/*}}}*/ +void Riftfront::GetNodesSidList(int* sidlist){/*{{{*/ + + _assert_(sidlist); + _assert_(nodes); + + for(int i=0;iSid(); +} +/*}}}*/ +int Riftfront::GetNumberOfNodes(void){/*{{{*/ + + return NUMVERTICES; +} +/*}}}*/ +bool Riftfront::InAnalysis(int in_analysis_type){/*{{{*/ + if (in_analysis_type==this->analysis_type) return true; + else return false; } /*}}}*/ bool Riftfront::IsPenalty(void){/*{{{*/ return true; -} -/*}}}*/ -void Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ - } /*}}}*/ @@ -306,38 +324,20 @@ } /*}}}*/ -void Riftfront::CreateKMatrix(Matrix* Kff, Matrix* Kfs){/*{{{*/ - /*do nothing: */ - return; -} -/*}}}*/ -void Riftfront::CreatePVector(Vector* pf){/*{{{*/ - /*do nothing: */ - return; -} -/*}}}*/ -void Riftfront::GetNodesSidList(int* sidlist){/*{{{*/ - - _assert_(sidlist); - _assert_(nodes); - - for(int i=0;iSid(); -} -/*}}}*/ -void Riftfront::GetNodesLidList(int* lidlist){/*{{{*/ - - _assert_(lidlist); - _assert_(nodes); - - for(int i=0;iLid(); -} -/*}}}*/ -int Riftfront::GetNumberOfNodes(void){/*{{{*/ - - return NUMVERTICES; -} -/*}}}*/ -bool Riftfront::InAnalysis(int in_analysis_type){/*{{{*/ - if (in_analysis_type==this->analysis_type) return true; - else return false; +void Riftfront::ResetHooks(){/*{{{*/ + + this->nodes=NULL; + this->elements=NULL; + this->matpar=NULL; + this->parameters=NULL; + + /*Get Element type*/ + this->hnodes->reset(); + this->helements->reset(); + this->hmatpar->reset(); + +} +/*}}}*/ +void Riftfront::SetCurrentConfiguration(Elements* elementsin,Loads* loadsin,Nodes* nodesin,Vertices* verticesin,Materials* materialsin,Parameters* parametersin){/*{{{*/ + } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 18925) +++ /issm/trunk-jpl/src/c/classes/Loads/Riftfront.h (revision 18926) @@ -57,42 +57,42 @@ /*}}}*/ /*Object virtual functions definitions:{{{ */ - void Echo(); - void DeepEcho(); - int Id(); - int ObjectEnum(); - Object* copy(); + Object* copy(); + void DeepEcho(); + void Echo(); + int Id(); + int ObjectEnum(); /*}}}*/ /*Update virtual functions resolution: {{{*/ - void InputUpdateFromVector(IssmDouble* vector, int name, int type); - void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows,int ncols, int name, int type){_error_("Not implemented yet!");} - void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");} void InputUpdateFromConstant(IssmDouble constant, int name); void InputUpdateFromConstant(int constant, int name){_error_("Not implemented yet!");} void InputUpdateFromConstant(bool constant, int name); void InputUpdateFromIoModel(int index, IoModel* iomodel){_error_("not implemented yet");}; + void InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows,int ncols, int name, int type){_error_("Not implemented yet!");} + void InputUpdateFromVector(IssmDouble* vector, int name, int type); + void InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){_error_("Not implemented yet!");} /*}}}*/ /*Load virtual functions definitions: {{{*/ void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); - void ResetHooks(); - void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); + void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; void CreateKMatrix(Matrix* Kff, Matrix* Kfs); void CreatePVector(Vector* pf); - void CreateJacobianMatrix(Matrix* Jff){_error_("Not implemented yet");}; + void GetNodesLidList(int* lidlist); void GetNodesSidList(int* sidlist); - void GetNodesLidList(int* lidlist); int GetNumberOfNodes(void); + bool InAnalysis(int analysis_type); bool IsPenalty(void); void PenaltyCreateJacobianMatrix(Matrix* Jff,IssmDouble kmax){_error_("Not implemented yet");}; void PenaltyCreateKMatrix(Matrix* Kff, Matrix* kfs, IssmDouble kmax); void PenaltyCreatePVector(Vector* pf, IssmDouble kmax); + void ResetHooks(); + void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); - bool InAnalysis(int analysis_type); /*}}}*/ /*Riftfront specific routines: {{{*/ + int Constrain(int* punstable); + void FreezeConstraints(void); + bool IsFrozen(void); ElementMatrix* PenaltyCreateKMatrixStressbalanceHoriz(IssmDouble kmax); ElementVector* PenaltyCreatePVectorStressbalanceHoriz(IssmDouble kmax); - int Constrain(int* punstable); - void FreezeConstraints(void); - bool IsFrozen(void); /*}}}*/ };