Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17000) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17001) @@ -61,6 +61,7 @@ void GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss); void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element); /*FS*/ + ElementVector* CreateDVectorFS(Element* element); ElementMatrix* CreateJacobianMatrixFS(Element* element); ElementMatrix* CreateKMatrixFS(Element* element); ElementMatrix* CreateKMatrixFSViscous(Element* element); Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17000) +++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17001) @@ -812,8 +812,17 @@ /*Finite Element Analysis*/ ElementVector* StressbalanceAnalysis::CreateDVector(Element* element){/*{{{*/ - /*Default, return NULL*/ + + int approximation; + element->GetInputValue(&approximation,ApproximationEnum); + switch(approximation){ + case FSApproximationEnum: + return CreateDVectorFS(element); + default: + return NULL; //no need for doftypes outside of FS approximation + } return NULL; + }/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ @@ -2380,6 +2389,35 @@ }/*}}}*/ /*FS*/ +ElementVector* StressbalanceAnalysis::CreateDVectorFS(Element* element){/*{{{*/ + + int meshtype,dim; + + /*Get problem dimension*/ + element->FindParam(&meshtype,MeshTypeEnum); + switch(meshtype){ + case Mesh2DverticalEnum: dim = 2; break; + case Mesh3DEnum: dim = 3; break; + default: _error_("mesh "<GetNumberOfNodesVelocity(); + int pnumnodes = element->GetNumberOfNodesPressure(); + + /*Initialize output vector*/ + ElementVector* de = element->NewElementVector(FSvelocityEnum); + + for(int i=0;ivalues[i*dim+j]=VelocityEnum; + } + for(int i=0;ivalues[vnumnodes*dim+i]=PressureEnum; + } + + return de; + +}/*}}}*/ ElementMatrix* StressbalanceAnalysis::CreateJacobianMatrixFS(Element* element){/*{{{*/ /*Intermediaries */ Index: ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 17000) +++ ../trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 17001) @@ -101,9 +101,13 @@ } /*Create dof vector for stiffness matrix preconditioning*/ - for (i=0;ielements->Size();i++){ - element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); - element->CreateDVector(df); + if(pdf){ + for(i=0;ielements->Size();i++){ + element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + ElementVector* de=analysis->CreateDVector(element); + if(de) de->InsertIntoGlobal(df); + delete de; + } } /*Assemble matrices and vector*/ Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17000) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 17001) @@ -103,7 +103,6 @@ virtual void Configure(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,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; - virtual void CreateDVector(Vector* df)=0; virtual void ElementSizes(IssmDouble* phx,IssmDouble* phy,IssmDouble* phz)=0; virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0; virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0; Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17000) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17001) @@ -190,12 +190,6 @@ return sqrt(2*this->GetArea()); } /*}}}*/ -/*FUNCTION Tria::CreateDVector {{{*/ -void Tria::CreateDVector(Vector* df){ - - /*Nothing done yet*/ -} -/*}}}*/ /*FUNCTION Tria::ComputeBasalStress {{{*/ void Tria::ComputeBasalStress(Vector* eps){ _error_("Not Implemented yet"); Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17000) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 17001) @@ -66,7 +66,6 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); - void CreateDVector(Vector* df); void Delta18oParameterization(void); void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17000) +++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17001) @@ -378,32 +378,6 @@ this->inputs->Configure(parameters); } /*}}}*/ -/*FUNCTION Penta::CreateDVector {{{*/ -void Penta::CreateDVector(Vector* df){ - - /*retrieve parameters: */ - ElementVector* De=NULL; - int analysis_type; - parameters->FindParam(&analysis_type,AnalysisTypeEnum); - - /*Checks in debugging*/ - _assert_(this->nodes && this->material && this->matpar && this->verticalneighbors && this->parameters && this->inputs); - - switch(analysis_type){ - #ifdef _HAVE_STRESSBALANCE_ - case StressbalanceAnalysisEnum: - De=CreateDVectorStressbalanceHoriz(); - break; - #endif - } - - /*Add to global Vector*/ - if(De){ - De->InsertIntoGlobal(df); - delete De; - } -} -/*}}}*/ /*FUNCTION Penta::DeepEcho{{{*/ void Penta::DeepEcho(void){ @@ -4475,52 +4449,6 @@ /*}}}*/ #endif -#ifdef _HAVE_STRESSBALANCE_ -/*FUNCTION Penta::CreateDVectorStressbalanceHoriz {{{*/ -ElementVector* Penta::CreateDVectorStressbalanceHoriz(void){ - - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - - switch(approximation){ - case FSApproximationEnum: - return CreateDVectorStressbalanceFS(); - default: - return NULL; //no need for doftypes outside of FS approximation - } -} -/*}}}*/ -/*FUNCTION Penta::CreateDVectorStressbalanceFS{{{*/ -ElementVector* Penta::CreateDVectorStressbalanceFS(void){ - - /*output: */ - ElementVector* De=NULL; - - /*Initialize Element vector and return if necessary*/ - int approximation; - inputs->GetInputValue(&approximation,ApproximationEnum); - if(approximation!=FSApproximationEnum) return NULL; - - /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocity(); - int pnumnodes = this->NumberofNodesPressure(); - - De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum); - - for(int i=0;ivalues[i*3+0]=VelocityEnum; - De->values[i*3+1]=VelocityEnum; - De->values[i*3+2]=VelocityEnum; - } - for(int i=0;ivalues[vnumnodes*3+i]=PressureEnum; - } - - return De; -} -/*}}}*/ -#endif - #ifdef _HAVE_HYDROLOGY_ /*FUNCTION Penta::CreateEPLDomainMassMatrix {{{*/ ElementMatrix* Penta::CreateEPLDomainMassMatrix(void){ Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17000) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 17001) @@ -70,7 +70,6 @@ int FiniteElement(void); void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters); void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum); - void CreateDVector(Vector* df); void Delta18oParameterization(void); void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure); Element* GetBasalElement(void); @@ -251,11 +250,6 @@ Tria* SpawnTria(int location); IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa); - #ifdef _HAVE_STRESSBALANCE_ - ElementVector* CreateDVectorStressbalanceHoriz(void); - ElementVector* CreateDVectorStressbalanceFS(void); - #endif - #ifdef _HAVE_HYDROLOGY_ ElementMatrix* CreateEPLDomainMassMatrix(void); void GetHydrologyDCInefficientHmax(IssmDouble* ph_max, Node* innode); Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17000) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 17001) @@ -68,7 +68,6 @@ void Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; void SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters){_error_("not implemented yet");}; void SetwiseNodeConnectivity(int* d_nz,int* o_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){_error_("not implemented yet");}; - void CreateDVector(Vector* df){_error_("not implemented yet");}; void Delta18oParameterization(void){_error_("not implemented yet");}; void ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");};