Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17000) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp (revision 17001) @@ -813,6 +813,15 @@ /*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){/*{{{*/ @@ -2381,4 +2390,33 @@ /*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){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17000) +++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h (revision 17001) @@ -62,4 +62,5 @@ void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element); /*FS*/ + ElementVector* CreateDVectorFS(Element* element); ElementMatrix* CreateJacobianMatrixFS(Element* element); ElementMatrix* CreateKMatrixFS(Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17000) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 17001) @@ -104,5 +104,4 @@ 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; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17000) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 17001) @@ -377,30 +377,4 @@ /*get inputs configured too: */ 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; - } } /*}}}*/ @@ -4476,50 +4450,4 @@ #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 {{{*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17000) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 17001) @@ -71,5 +71,4 @@ 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); @@ -252,9 +251,4 @@ 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); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17000) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 17001) @@ -69,5 +69,4 @@ 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");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17000) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 17001) @@ -189,10 +189,4 @@ return sqrt(2*this->GetArea()); -} -/*}}}*/ -/*FUNCTION Tria::CreateDVector {{{*/ -void Tria::CreateDVector(Vector* df){ - - /*Nothing done yet*/ } /*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17000) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 17001) @@ -67,5 +67,4 @@ 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");}; Index: /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 17000) +++ /issm/trunk-jpl/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 17001) @@ -102,7 +102,11 @@ /*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; + } }