Index: /issm/trunk/src/c/Container/Constraints.cpp =================================================================== --- /issm/trunk/src/c/Container/Constraints.cpp (revision 5771) +++ /issm/trunk/src/c/Container/Constraints.cpp (revision 5772) @@ -63,45 +63,2 @@ } /*}}}*/ -/*FUNCTION Constraints::SetupSpcs{{{1*/ -void Constraints::SetupSpcs(Nodes* nodes,Vec yg,int analysis_type){ - - int i; - - Node* node=NULL; - int nodeid; - int dof; - double value; - - for(i=0;iSize();i++){ - - Object* object=(Object*)this->GetObjectByOffset(i); - - /*Check this is a single point constraint (spc): */ - if(object->Enum()==SpcEnum){ - - Spc* spc=(Spc*)object; - - if(spc->InAnalysis(analysis_type)){ - - /*Ok, this object is a constraint. Get the nodeid from the node it applies to: */ - nodeid=spc->GetNodeId(); - dof=spc->GetDof(); - value=spc->GetValue(); - - /*Now, chase through nodes and find the corect node: */ - node=(Node*)nodes->GetObjectById(NULL,nodeid); - - /*Apply constraint: */ - if(node){ //in case the spc is dealing with a node on another cpu - node->ApplyConstraint(yg,dof,value); - } - } - - } - } - - /*Assemble yg: */ - VecAssemblyBegin(yg); - VecAssemblyEnd(yg); -} -/*}}}*/ Index: /issm/trunk/src/c/Container/Constraints.h =================================================================== --- /issm/trunk/src/c/Container/Constraints.h (revision 5771) +++ /issm/trunk/src/c/Container/Constraints.h (revision 5772) @@ -28,5 +28,4 @@ /*numerics: {{{1*/ int NumberOfConstraints(void); - void SetupSpcs(Nodes* nodes,Vec yg,int analysis_type); /*}}}*/ Index: /issm/trunk/src/c/Container/Nodes.cpp =================================================================== --- /issm/trunk/src/c/Container/Nodes.cpp (revision 5771) +++ /issm/trunk/src/c/Container/Nodes.cpp (revision 5772) @@ -44,5 +44,5 @@ /*Numerics*/ /*FUNCTION Nodes::DistributeDofs{{{1*/ -void Nodes::DistributeDofs(int analysis_type){ +void Nodes::DistributeDofs(int analysis_type,int setenum){ extern int num_procs; @@ -58,4 +58,7 @@ int numnodes=0; + /*some check: */ + if ((setenum!=GsetEnum) && (setenum!=FsetEnum) && (setenum!=SsetEnum))ISSMERROR("%s%s%s"," dof distribution for set of enum type ",EnumToString(setenum)," not supported yet!"); + /*Go through objects, and distribute dofs locally, from 0 to numberofdofs: */ for (i=0;iSize();i++){ @@ -64,9 +67,7 @@ /*Check that this node corresponds to our analysis currently being carried out: */ if (node->InAnalysis(analysis_type)){ - node->DistributeDofs(&dofcount); - } - } - - + node->DistributeDofs(&dofcount,setenum); + } + } /*Ok, now every object has distributed dofs, but locally, and with a dof count starting from @@ -96,5 +97,5 @@ Node* node=(Node*)this->GetObjectByOffset(i); if (node->InAnalysis(analysis_type)){ - node->OffsetDofs(dofcount); + node->OffsetDofs(dofcount,setenum); } @@ -104,5 +105,5 @@ * object that is not a clone, tell them to show their dofs, so that later on, they can get picked * up by their clones: */ - maxdofspernode=this->MaxNumDofs(analysis_type); + maxdofspernode=this->MaxNumDofs(analysis_type,setenum); numnodes=this->NumberOfNodes(analysis_type); @@ -114,5 +115,5 @@ Node* node=(Node*)this->GetObjectByOffset(i); if (node->InAnalysis(analysis_type)){ - node->ShowTrueDofs(truedofs,maxdofspernode);//give maxdofspernode, column size, so that nodes can index into truedofs + node->ShowTrueDofs(truedofs,maxdofspernode,setenum);//give maxdofspernode, column size, so that nodes can index into truedofs } } @@ -125,5 +126,5 @@ Node* node=(Node*)this->GetObjectByOffset(i); if (node->InAnalysis(analysis_type)){ - node->UpdateCloneDofs(alltruedofs,maxdofspernode); //give maxdofspernode, column size, so that nodes can index into alltruedofs + node->UpdateCloneDofs(alltruedofs,maxdofspernode,setenum); //give maxdofspernode, column size, so that nodes can index into alltruedofs } } @@ -215,5 +216,5 @@ /*}}}*/ /*FUNCTION Nodes::NumberOfDofs{{{1*/ -int Nodes::NumberOfDofs(int analysis_type){ +int Nodes::NumberOfDofs(int analysis_type,int setenum){ int i; @@ -233,5 +234,5 @@ if (!node->IsClone()){ - numdofs+=node->GetNumberOfDofs(); + numdofs+=node->GetNumberOfDofs(NoneApproximationEnum,setenum); } @@ -241,5 +242,4 @@ /*Gather from all cpus: */ MPI_Allreduce ( (void*)&numdofs,(void*)&allnumdofs,1,MPI_INT,MPI_SUM,MPI_COMM_WORLD); - return allnumdofs; } @@ -335,5 +335,5 @@ /*}}}*/ /*FUNCTION Nodes::MaxNumDofs{{{1*/ -int Nodes::MaxNumDofs(int analysis_type){ +int Nodes::MaxNumDofs(int analysis_type,int setenum){ int i; @@ -350,5 +350,5 @@ if (node->InAnalysis(analysis_type)){ - numdofs=node->GetNumberOfDofs(); + numdofs=node->GetNumberOfDofs(NoneApproximationEnum,setenum); if (numdofs>max)max=numdofs; } Index: /issm/trunk/src/c/Container/Nodes.h =================================================================== --- /issm/trunk/src/c/Container/Nodes.h (revision 5771) +++ /issm/trunk/src/c/Container/Nodes.h (revision 5772) @@ -19,12 +19,12 @@ /*}}}*/ /*numerics: {{{1*/ - void DistributeDofs(int analysis_type); + void DistributeDofs(int analysis_type,int SETENUM); void FlagClones(int analysis_type); void FlagNodeSets(Vec pv_g, Vec pv_f, Vec pv_s,int analysis_type); - int NumberOfDofs(int analysis_type); + int NumberOfDofs(int analysis_type,int setenum); int NumberOfNodes(int analysis_type); int NumberOfNodes(void); void Ranks(int* ranks,int analysis_type); - int MaxNumDofs(int analysis_type); + int MaxNumDofs(int analysis_type,int setenum); /*}}}*/ Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 5771) +++ /issm/trunk/src/c/Makefile.am (revision 5772) @@ -172,4 +172,8 @@ ./objects/Loads/Numericalflux.cpp\ ./objects/Loads/Numericalflux.h\ + ./objects/Numerics/ElementMatrix.h\ + ./objects/Numerics/ElementMatrix.cpp\ + ./objects/Numerics/ElementVector.h\ + ./objects/Numerics/ElementVector.cpp\ ./objects/Params/Param.h\ ./objects/Params/BoolParam.cpp\ @@ -448,4 +452,6 @@ ./modules/CostFunctionx/CostFunctionx.h\ ./modules/CostFunctionx/CostFunctionx.cpp\ + ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\ + ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\ ./modules/Orthx/Orthx.h\ ./modules/Orthx/Orthx.cpp\ @@ -524,4 +530,6 @@ ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h\ ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\ + ./modules/Reduceloadx/Reduceloadx.h\ + ./modules/Reduceloadx/Reduceloadx.cpp\ ./modules/NodeConnectivityx/NodeConnectivityx.cpp\ ./modules/NodeConnectivityx/NodeConnectivityx.h\ @@ -723,4 +731,8 @@ ./objects/Loads/Numericalflux.cpp\ ./objects/Loads/Numericalflux.h\ + ./objects/Numerics/ElementMatrix.h\ + ./objects/Numerics/ElementMatrix.cpp\ + ./objects/Numerics/ElementVector.h\ + ./objects/Numerics/ElementVector.cpp\ ./objects/Params/Param.h\ ./objects/Params/BoolParam.cpp\ @@ -996,4 +1008,6 @@ ./modules/CostFunctionx/CostFunctionx.h\ ./modules/CostFunctionx/CostFunctionx.cpp\ + ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h\ + ./modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp\ ./modules/Orthx/Orthx.h\ ./modules/Orthx/Orthx.cpp\ @@ -1070,4 +1084,6 @@ ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.h\ ./modules/Reduceloadfromgtofx/Reduceloadfromgtofx.cpp\ + ./modules/Reduceloadx/Reduceloadx.h\ + ./modules/Reduceloadx/Reduceloadx.cpp\ ./modules/MassFluxx/MassFluxx.cpp\ ./modules/MassFluxx/MassFluxx.h\ @@ -1128,4 +1144,6 @@ ./solvers/solver_adjoint_linear.cpp\ ./solvers/solver_diagnostic_nonlinear.cpp\ + ./solvers/solver_diagnostic_nonlinear_kgg.cpp\ + ./solvers/solver_diagnostic_nonlinear_kff.cpp\ ./solvers/solver_thermal_nonlinear.cpp\ ./modules/Bamgx/Bamgx.cpp\ Index: /issm/trunk/src/c/modules/BuildNodeSetsx/BuildNodeSetsx.cpp =================================================================== --- /issm/trunk/src/c/modules/BuildNodeSetsx/BuildNodeSetsx.cpp (revision 5771) +++ /issm/trunk/src/c/modules/BuildNodeSetsx/BuildNodeSetsx.cpp (revision 5772) @@ -29,5 +29,5 @@ /*Get gsize; */ - gsize=nodes->NumberOfDofs(analysis_type); + gsize=nodes->NumberOfDofs(analysis_type,GsetEnum); if(gsize){ Index: /issm/trunk/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp =================================================================== --- /issm/trunk/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp (revision 5772) +++ /issm/trunk/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.cpp (revision 5772) @@ -0,0 +1,43 @@ +/*!\file CreateNodalConstraintsx + * \brief: establish degrees of freedom for all nodes, and return partitioning vector. Do only once. + */ + +#include "./CreateNodalConstraintsx.h" + +#include "../../shared/shared.h" +#include "../../include/include.h" +#include "../../toolkits/toolkits.h" +#include "../../EnumDefinitions/EnumDefinitions.h" + +void CreateNodalConstraintsx( Vec* pys, Nodes* nodes,int analysis_type){ + + int i; + + /*intermediary: */ + int numberofdofs; + + /*output: */ + Vec ys=NULL; + + /*figure out how many dofs we have: */ + numberofdofs=nodes->NumberOfDofs(analysis_type,SsetEnum); + + /*allocate:*/ + ys=NewVec(numberofdofs); + + /*go through all nodes, and for the ones corresponding to this analysis_type, fill the + * constraints vector with the constraint values: */ + for (i=0;iSize();i++){ + Node* node=(Node*)nodes->GetObjectByOffset(i); + if (node->InAnalysis(analysis_type)){ + node->CreateNodalConstraints(ys); + } + } + + /*Assemble: */ + VecAssemblyBegin(ys); + VecAssemblyEnd(ys); + + /*Assign output pointers: */ + *pys=ys; +} Index: /issm/trunk/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h =================================================================== --- /issm/trunk/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h (revision 5772) +++ /issm/trunk/src/c/modules/CreateNodalConstraintsx/CreateNodalConstraintsx.h (revision 5772) @@ -0,0 +1,14 @@ +/*!\file: CreateNodalConstraintsx.h + */ + +#ifndef _CREATENODALCONSTRAINTSX_H +#define _CREATENODALCONSTRAINTSX_H + +#include "../../Container/Container.h" +#include "../../objects/objects.h" + +/* local prototypes: */ +void CreateNodalConstraintsx( Vec* pys, Nodes* nodes,int analysis_type); + +#endif /* _CREATENODALCONSTRAINTSX_H */ + Index: /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp =================================================================== --- /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp (revision 5771) +++ /issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.cpp (revision 5772) @@ -25,5 +25,5 @@ /*Get size of vector: */ - gsize=nodes->NumberOfDofs(configuration_type); + gsize=nodes->NumberOfDofs(configuration_type,GsetEnum); if (gsize==0) ISSMERROR("Allocating a Vec of size 0 as gsize=0 for configuration: %s",EnumToString(configuration_type)); Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 5771) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 5772) @@ -66,4 +66,5 @@ parameters->AddObject(new StringParam(SolverStringEnum,iomodel->solverstring)); parameters->AddObject(new IntParam(NumberOfElementsEnum,iomodel->numberofelements)); + parameters->AddObject(new BoolParam(KffEnum,iomodel->kff)); /*Deal with more complex parameters*/ Index: /issm/trunk/src/c/modules/NodesDofx/NodesDofx.cpp =================================================================== --- /issm/trunk/src/c/modules/NodesDofx/NodesDofx.cpp (revision 5771) +++ /issm/trunk/src/c/modules/NodesDofx/NodesDofx.cpp (revision 5772) @@ -26,6 +26,7 @@ *a node has already been distributed dofs on one cpu, all other cpus with the same node cannot distribute it *anymore. Use clone field to be sure of that: */ - nodes->DistributeDofs(analysis_type); - + nodes->DistributeDofs(analysis_type,GsetEnum); + nodes->DistributeDofs(analysis_type,FsetEnum); + nodes->DistributeDofs(analysis_type,SsetEnum); } Index: /issm/trunk/src/c/modules/Reduceloadx/Reduceloadx.cpp =================================================================== --- /issm/trunk/src/c/modules/Reduceloadx/Reduceloadx.cpp (revision 5772) +++ /issm/trunk/src/c/modules/Reduceloadx/Reduceloadx.cpp (revision 5772) @@ -0,0 +1,49 @@ +/*!\file Reduceloadx + * \brief reduce loads (wring out boundary conditions) + */ + +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./Reduceloadx.h" + +void Reduceloadx( Vec pf, Mat Kfs, Vec y_s, bool flag_ys0){ + + /*intermediary*/ + Vec y_s0=NULL; + Vec Kfsy_s=NULL; + int Kfsm,Kfsn; + PetscScalar a; + + if(pf){ + + /*pf = pf - Kfs * y_s;*/ + MatGetLocalSize(Kfs,&Kfsm,&Kfsn); + Kfsy_s=NewVecFromLocalSize(Kfsm); + if (flag_ys0){ + + /*Create y_s0, full of 0: */ + VecDuplicate(y_s,&y_s0); + VecSet(y_s0,0.0); + VecAssemblyBegin(y_s0); + VecAssemblyEnd(y_s0); + + MatMultPatch(Kfs,y_s0,Kfsy_s); + } + else{ + MatMultPatch(Kfs,y_s,Kfsy_s); + } + + a=-1; + VecAXPY(pf,a,Kfsy_s); + } + + + /*Free ressources and return*/ + VecFree(&y_s0); + VecFree(&Kfsy_s); + +} Index: /issm/trunk/src/c/modules/Reduceloadx/Reduceloadx.h =================================================================== --- /issm/trunk/src/c/modules/Reduceloadx/Reduceloadx.h (revision 5772) +++ /issm/trunk/src/c/modules/Reduceloadx/Reduceloadx.h (revision 5772) @@ -0,0 +1,13 @@ +/*!\file: Reduceloadx.h + * \brief reduce loads (wring out boundary conditions) + */ + +#ifndef _REDUCELOADX_H +#define _REDUCELOADX_H + +#include "../../objects/objects.h" + +/* local prototypes: */ +void Reduceloadx( Vec pf, Mat Kfs, Vec ys, bool flag_ys0=false); + +#endif /* _REDUCELOADX_H */ Index: /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp =================================================================== --- /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp (revision 5771) +++ /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.cpp (revision 5772) @@ -10,29 +10,39 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void SpcNodesx(Vec* pyg, Nodes* nodes,Constraints* constraints,int analysis_type){ +void SpcNodesx(Nodes* nodes,Constraints* constraints,int analysis_type){ int i; - int numberofdofs; - int gsize; + Node* node=NULL; + int nodeid; + int dof; + double value; - /*output: */ - Vec yg=NULL; + for(i=0;iSize();i++){ - /*First, recover number of dofs from nodes: */ - numberofdofs=nodes->NumberOfDofs(analysis_type); + Object* object=(Object*)constraints->GetObjectByOffset(i); - if(numberofdofs){ - - /*Allocate yg: */ - yg=NewVec(numberofdofs); + /*Check constraints is a single point constraint (spc): */ + if(object->Enum()==SpcEnum){ - /*Now, go through constraints, and update the nodes and the constraint vector at the same time: */ - constraints->SetupSpcs(nodes,yg,analysis_type); + Spc* spc=(Spc*)object; - /*Specify numentries: */ - VecGetSize(yg,&gsize); + if(spc->InAnalysis(analysis_type)){ + + /*Ok, constraints object is a constraint. Get the nodeid from the node it applies to: */ + nodeid=spc->GetNodeId(); + dof=spc->GetDof(); + value=spc->GetValue(); + + /*Now, chase through nodes and find the corect node: */ + node=(Node*)nodes->GetObjectById(NULL,nodeid); + + /*Apply constraint: */ + if(node){ //in case the spc is dealing with a node on another cpu + node->ApplyConstraint(dof,value); + } + } + + } } - /*Assign output pointers: */ - *pyg=yg; } Index: /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h =================================================================== --- /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h (revision 5771) +++ /issm/trunk/src/c/modules/SpcNodesx/SpcNodesx.h (revision 5772) @@ -11,5 +11,5 @@ /* local prototypes: */ -void SpcNodesx(Vec* pyg, Nodes* nodesin,Constraints* constraints,int analysis_type); +void SpcNodesx(Nodes* nodesin,Constraints* constraints,int analysis_type); #endif /* _SPCNODESX_H */ Index: /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp =================================================================== --- /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 5771) +++ /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.cpp (revision 5772) @@ -9,16 +9,21 @@ #include "../../EnumDefinitions/EnumDefinitions.h" -void SystemMatricesx(Mat* pKgg, Vec* ppg,double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){ +void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters,bool kflag,bool pflag,bool penalty_kflag,bool penalty_pflag){ /*intermediary: */ - int i,j,gsize; + int i,j; + int gsize,fsize,ssize; int connectivity, numberofdofspernode; int analysis_type, configuration_type; Element *element = NULL; Load *load = NULL; + bool buildkff=false; /*output: */ Mat Kgg = NULL; + Mat Kff = NULL; + Mat Kfs = NULL; Vec pg = NULL; + Vec pf = NULL; double kmax = 0; @@ -31,8 +36,11 @@ parameters->FindParam(&configuration_type,ConfigurationTypeEnum); parameters->FindParam(&connectivity,ConnectivityEnum); + parameters->FindParam(&buildkff,KffEnum); /*Get size of matrix: */ - gsize=nodes->NumberOfDofs(configuration_type); - numberofdofspernode=nodes->MaxNumDofs(configuration_type); + gsize=nodes->NumberOfDofs(configuration_type,GsetEnum); + fsize=nodes->NumberOfDofs(configuration_type,FsetEnum); + ssize=nodes->NumberOfDofs(configuration_type,SsetEnum); + numberofdofspernode=nodes->MaxNumDofs(configuration_type,GsetEnum); /*Checks in debugging mode {{{1*/ @@ -46,10 +54,14 @@ if(kflag){ - Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode); + if(!buildkff)Kgg=NewMat(gsize,gsize,NULL,&connectivity,&numberofdofspernode); + else{ + Kff=NewMat(fsize,fsize,NULL,&connectivity,&numberofdofspernode); + Kfs=NewMat(fsize,ssize,NULL,&connectivity,&numberofdofspernode); + } /*Fill stiffness matrix from elements: */ for (i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - element->CreateKMatrix(Kgg); + element->CreateKMatrix(Kgg,Kff,Kfs); } @@ -57,11 +69,22 @@ for (i=0;iSize();i++){ load=(Load*)loads->GetObjectByOffset(i); - if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kgg); + if (load->InAnalysis(configuration_type)) load->CreateKMatrix(Kgg,Kff,Kfs); } /*Assemble matrix and compress matrix to save memory: */ - MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY); - MatCompress(Kgg); + if(!buildkff){ + MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY); + MatCompress(Kgg); + } + else{ + MatAssemblyBegin(Kff,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kff,MAT_FINAL_ASSEMBLY); + MatCompress(Kff); + + MatAssemblyBegin(Kfs,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kfs,MAT_FINAL_ASSEMBLY); + MatCompress(Kfs); + } } if(pflag){ @@ -72,5 +95,5 @@ for (i=0;iSize();i++){ element=(Element*)elements->GetObjectByOffset(i); - element->CreatePVector(pg); + element->CreatePVector(pg,pf); } @@ -78,27 +101,46 @@ for (i=0;iSize();i++){ load=(Load*)loads->GetObjectByOffset(i); - if (load->InAnalysis(configuration_type)) load->CreatePVector(pg); + if (load->InAnalysis(configuration_type)) load->CreatePVector(pg,pf); } - VecAssemblyBegin(pg); - VecAssemblyEnd(pg); + if(!buildkff){ + VecAssemblyBegin(pg); + VecAssemblyEnd(pg); + } + else{ + VecAssemblyBegin(pf); + VecAssemblyEnd(pf); + } } + /*Now, deal with penalties*/ if(penalty_kflag){ /*Now, figure out maximum value of K_gg, so that we can penalize it correctly: */ - MatNorm(Kgg,NORM_INFINITY,&kmax); + if(!buildkff)MatNorm(Kgg,NORM_INFINITY,&kmax); + else MatNorm(Kff,NORM_INFINITY,&kmax); /*Fill stiffness matrix from loads: */ for (i=0;iSize();i++){ load=(Load*)loads->GetObjectByOffset(i); - if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kgg,kmax); + if (load->InAnalysis(configuration_type)) load->PenaltyCreateKMatrix(Kgg,Kff,Kfs,kmax); } /*Assemble matrix and compress matrix to save memory: */ - MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY); - MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY); - MatCompress(Kgg); + if(!buildkff){ + MatAssemblyBegin(Kgg,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kgg,MAT_FINAL_ASSEMBLY); + MatCompress(Kgg); + } + else{ + MatAssemblyBegin(Kff,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kff,MAT_FINAL_ASSEMBLY); + MatCompress(Kff); + + MatAssemblyBegin(Kfs,MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(Kfs,MAT_FINAL_ASSEMBLY); + MatCompress(Kfs); + } } if(penalty_pflag){ @@ -107,14 +149,23 @@ for (i=0;iSize();i++){ load=(Load*)loads->GetObjectByOffset(i); - if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pg,kmax); + if (load->InAnalysis(configuration_type)) load->PenaltyCreatePVector(pg,pf,kmax); } - VecAssemblyBegin(pg); - VecAssemblyEnd(pg); + if(!buildkff){ + VecAssemblyBegin(pg); + VecAssemblyEnd(pg); + } + else{ + VecAssemblyBegin(pf); + VecAssemblyEnd(pf); + } } /*Assign output pointers: */ - *pKgg=Kgg; - *ppg=pg; + if(pKgg) *pKgg=Kgg; + if(pKff) *pKff=Kff; + if(pKfs) *pKfs=Kfs; + if(ppg) *ppg=pg; + if(ppf) *ppf=pf; if(pkmax) *pkmax=kmax; } Index: /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h =================================================================== --- /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h (revision 5771) +++ /issm/trunk/src/c/modules/SystemMatricesx/SystemMatricesx.h (revision 5772) @@ -10,5 +10,5 @@ /* local prototypes: */ -void SystemMatricesx(Mat* pKgg, Vec* ppg,double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, +void SystemMatricesx(Mat* pKgg, Mat* pKff, Mat* pKfs, Vec* ppg, Vec* ppf, double* pkmax,Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads,Materials* materials, Parameters* parameters, bool kflag=true,bool pflag=true,bool penalty_kflag=true,bool penalty_pflag=true); Index: /issm/trunk/src/c/modules/modules.h =================================================================== --- /issm/trunk/src/c/modules/modules.h (revision 5771) +++ /issm/trunk/src/c/modules/modules.h (revision 5772) @@ -20,4 +20,5 @@ #include "./ContourToNodesx/ContourToNodesx.h" #include "./CostFunctionx/CostFunctionx.h" +#include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h" #include "./DakotaResponsesx/DakotaResponsesx.h" #include "./ElementConnectivityx/ElementConnectivityx.h" @@ -72,4 +73,5 @@ #include "./Qmux/Qmux.h" #include "./Reduceloadfromgtofx/Reduceloadfromgtofx.h" +#include "./Reduceloadx/Reduceloadx.h" #include "./Reducematrixfromgtofx/Reducematrixfromgtofx.h" #include "./Reducevectorgtosx/Reducevectorgtosx.h" Index: /issm/trunk/src/c/objects/DofIndexing.cpp =================================================================== --- /issm/trunk/src/c/objects/DofIndexing.cpp (revision 5771) +++ /issm/trunk/src/c/objects/DofIndexing.cpp (revision 5772) @@ -21,16 +21,21 @@ DofIndexing::DofIndexing(){ - int i; - this->numberofdofs=UNDEF; + this->gsize=UNDEF; + this->fsize=UNDEF; + this->ssize=UNDEF; this->clone=0; this->f_set=NULL; this->s_set=NULL; - this->doflist=NULL; + this->svalues=NULL; this->doftype=NULL; -} -/*}}}*/ -/*FUNCTION DofIndexing::DofIndexing(int numberofdofs){{{1*/ -DofIndexing::DofIndexing(int in_numberofdofs){ - this->Init(in_numberofdofs,NULL); + this->gdoflist=NULL; + this->fdoflist=NULL; + this->sdoflist=NULL; + +} +/*}}}*/ +/*FUNCTION DofIndexing::DofIndexing(int gsize){{{1*/ +DofIndexing::DofIndexing(int in_gsize){ + this->Init(in_gsize,NULL); } /*}}}*/ @@ -39,19 +44,37 @@ int i; - this->numberofdofs=in->numberofdofs; + this->gsize=in->gsize; + this->fsize=in->fsize; + this->ssize=in->ssize; + this->clone=in->clone; - this->f_set=(int*)xmalloc(this->numberofdofs*sizeof(int)); - this->s_set=(int*)xmalloc(this->numberofdofs*sizeof(int)); - this->doflist=(int*)xmalloc(this->numberofdofs*sizeof(int)); - if (in->doftype) this->doftype=(int*)xmalloc(this->numberofdofs*sizeof(int)); - else this->doftype=NULL; - - for(i=0;inumberofdofs;i++){ - this->f_set[i]=in->f_set[i]; - this->s_set[i]=in->s_set[i]; - this->doflist[i]=in->doflist[i]; - if (in->doftype) this->doftype[i]=in->doftype[i]; - } + if(this->gsize){ + this->f_set=(bool*)xmalloc(this->gsize*sizeof(bool)); + this->s_set=(bool*)xmalloc(this->gsize*sizeof(bool)); + this->svalues=(double*)xmalloc(this->gsize*sizeof(int)); + if(in->doftype)this->doftype=(int*)xmalloc(this->gsize*sizeof(int)); + this->gdoflist=(int*)xmalloc(this->gsize*sizeof(int)); + } + else{ + this->f_set=NULL; + this->s_set=NULL; + this->svalues=NULL; + this->doftype=NULL; + this->gdoflist=NULL; + } + if(this->fsize)this->fdoflist=(int*)xmalloc(this->fsize*sizeof(int)); else this->fdoflist=NULL; + if(this->ssize)this->sdoflist=(int*)xmalloc(this->ssize*sizeof(int)); else this->sdoflist=NULL; + + if(this->gsize){ + memcpy(this->f_set,in->f_set,this->gsize*sizeof(bool)); + memcpy(this->s_set,in->s_set,this->gsize*sizeof(bool)); + memcpy(this->svalues,in->svalues,this->gsize*sizeof(double)); + if(this->doftype)memcpy(this->doftype,in->doftype,this->gsize*sizeof(int)); + memcpy(this->gdoflist,in->gdoflist,this->gsize*sizeof(int)); + } + if(this->fsize)memcpy(this->fdoflist,in->fdoflist,this->fsize*sizeof(int)); + if(this->ssize)memcpy(this->sdoflist,in->sdoflist,this->ssize*sizeof(int)); + } /*}}}*/ @@ -61,29 +84,63 @@ xfree((void**)&f_set); xfree((void**)&s_set); - xfree((void**)&doflist); + xfree((void**)&svalues); xfree((void**)&doftype); + xfree((void**)&gdoflist); + xfree((void**)&fdoflist); + xfree((void**)&sdoflist); } /*}}}*/ /*FUNCTION DofIndexing::Init{{{1*/ -void DofIndexing::Init(int in_numberofdofs,int* in_doftype){ - - int i; - this->numberofdofs=in_numberofdofs; +void DofIndexing::Init(int in_gsize,int* in_doftype){ + + int i; + this->gsize=in_gsize; + this->clone=0; /*allocate: */ - this->f_set=(int*)xmalloc(this->numberofdofs*sizeof(int)); - this->s_set=(int*)xmalloc(this->numberofdofs*sizeof(int)); - this->doflist=(int*)xmalloc(this->numberofdofs*sizeof(int)); - if (in_doftype) this->doftype=(int*)xmalloc(this->numberofdofs*sizeof(int)); - - for (i=0;inumberofdofs;i++){ + if(this->gsize){ + this->f_set=(bool*)xmalloc(this->gsize*sizeof(bool)); + this->s_set=(bool*)xmalloc(this->gsize*sizeof(bool)); + this->svalues=(double*)xmalloc(this->gsize*sizeof(double)); + if(in_doftype)this->doftype=(int*)xmalloc(this->gsize*sizeof(int)); + this->gdoflist=(int*)xmalloc(this->gsize*sizeof(int)); + } + + for (i=0;igsize;i++){ /*assume dof is free, no constraints, no rigid body constraint: */ - this->f_set[i]=1; - this->s_set[i]=0; - this->doflist[i]=UNDEF; - if(in_doftype) this->doftype[i]=in_doftype[i]; - } + this->f_set[i]=true; + this->s_set[i]=false; + if(this->doftype)this->doftype[i]=in_doftype[i]; + this->svalues[i]=0; //0 constraint is the default value + } +} +/*}}}*/ +/*FUNCTION DofIndexing::InitSet{{{1*/ +void DofIndexing::InitSet(int setenum){ + + int i; + int size=0; + + /*go through sets, and figure out how many dofs belong to this set, except for g-set, + * which has already been initialized: */ + if(setenum==FsetEnum){ + size=0; + for(i=0;igsize;i++) if(f_set[i])size++; + this->fsize=size; + if(this->fsize)this->fdoflist=(int*)xmalloc(size*sizeof(int)); + else this->fdoflist=NULL; + } + else if(setenum==SsetEnum){ + size=0; + for(i=0;igsize;i++) if(s_set[i])size++; + this->ssize=size; + if(this->ssize)this->sdoflist=(int*)xmalloc(size*sizeof(int)); + else this->sdoflist=NULL; + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + + } /*}}}*/ @@ -96,5 +153,5 @@ printf("DofIndexing:\n"); - printf(" numberofdofs: %i\n",numberofdofs); + printf(" gsize: %i\n",gsize); printf(" clone: %i\n",clone); } @@ -106,16 +163,17 @@ printf("DofIndexing:\n"); - printf(" numberofdofs: %i\n",numberofdofs); + printf(" gsize: %i\n",gsize); + printf(" fsize: %i\n",fsize); + printf(" ssize: %i\n",ssize); printf(" clone: %i\n",clone); printf(" set membership: f,s sets \n"); - for(i=0;issize); + for(i=0;igsize;i++){ + if(this->s_set[i])printf(" %g |",svalues[i]); } printf("\n"); @@ -123,11 +181,28 @@ if(doftype){ printf(" doftype: |"); - for(i=0;igsize); + for(i=0;igsize;i++){ + printf(" %i |",gdoflist[i]); + } + printf("\n"); + + printf(" f_doflist (%i): |",this->fsize); + for(i=0;ifsize;i++){ + printf(" %i |",fdoflist[i]); + } + printf("\n"); + + printf(" s_doflist (%i): |",this->ssize); + for(i=0;issize;i++){ + printf(" %i |",sdoflist[i]); + } + printf("\n"); } /*}}}*/ @@ -145,22 +220,43 @@ memcpy(&enum_type,marshalled_dataset,sizeof(int)); marshalled_dataset+=sizeof(int); - memcpy(&numberofdofs,marshalled_dataset,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs); + /*easy part: */ + memcpy(&gsize,marshalled_dataset,sizeof(gsize));marshalled_dataset+=sizeof(gsize); + memcpy(&fsize,marshalled_dataset,sizeof(fsize));marshalled_dataset+=sizeof(fsize); + memcpy(&ssize,marshalled_dataset,sizeof(ssize));marshalled_dataset+=sizeof(ssize); + memcpy(&flagdoftype,marshalled_dataset,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype); memcpy(&clone,marshalled_dataset,sizeof(clone));marshalled_dataset+=sizeof(clone); /*Allocate: */ - this->f_set=(int*)xmalloc(this->numberofdofs*sizeof(int)); - this->s_set=(int*)xmalloc(this->numberofdofs*sizeof(int)); - this->doflist=(int*)xmalloc(this->numberofdofs*sizeof(int)); - - memcpy(f_set,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - memcpy(s_set,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - memcpy(doflist,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - - memcpy(&flagdoftype,marshalled_dataset,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype); - if(flagdoftype){ // there is a hook so demarshall it - this->doftype=(int*)xmalloc(this->numberofdofs*sizeof(int)); - memcpy(doftype,marshalled_dataset,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - } - else this->doftype=NULL; //There is no coupling, so doftype is NULL + if(this->gsize){ + this->f_set=(bool*)xmalloc(this->gsize*sizeof(bool)); + this->s_set=(bool*)xmalloc(this->gsize*sizeof(bool)); + this->svalues=(double*)xmalloc(this->ssize*sizeof(double)); + if(flagdoftype)this->doftype=(int*)xmalloc(this->gsize*sizeof(int)); + else this->doftype=NULL; + this->gdoflist=(int*)xmalloc(this->gsize*sizeof(int)); + } + else{ + this->f_set=NULL; + this->s_set=NULL; + this->svalues=NULL; + this->doftype=NULL; + this->gdoflist=NULL; + } + if(this->fsize)this->fdoflist=(int*)xmalloc(this->fsize*sizeof(int)); + else this->fdoflist=NULL; + if(this->ssize)this->sdoflist=(int*)xmalloc(this->ssize*sizeof(int)); + else this->sdoflist=NULL; + + /*Copy arrays: */ + if(this->gsize){ + memcpy(f_set,marshalled_dataset,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool); + memcpy(s_set,marshalled_dataset,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool); + memcpy(this->svalues,marshalled_dataset,this->gsize*sizeof(double));marshalled_dataset+=this->gsize*sizeof(double); + if(flagdoftype){ memcpy(doftype,marshalled_dataset,gsize*sizeof(int));marshalled_dataset+=gsize*sizeof(int); } + memcpy(this->gdoflist,marshalled_dataset,this->gsize*sizeof(int));marshalled_dataset+=this->gsize*sizeof(int); + } + + if(this->fsize){ memcpy(this->fdoflist,marshalled_dataset,this->fsize*sizeof(int));marshalled_dataset+=this->fsize*sizeof(int); } + if(this->ssize){ memcpy(this->sdoflist,marshalled_dataset,this->ssize*sizeof(int));marshalled_dataset+=this->ssize*sizeof(int); } /*return: */ @@ -174,9 +270,13 @@ char* marshalled_dataset=NULL; int enum_type=0; - int flagdoftype; //to indicate if there are some doftype or if NULL + bool flagdoftype; //to indicate if there are some doftype or if NULL /*recover marshalled_dataset: */ marshalled_dataset=*pmarshalled_dataset; + /*preliminary: */ + if(this->doftype)flagdoftype=true; + else flagdoftype=false; + /*get enum type of DofIndexing: */ enum_type=DofIndexingEnum; @@ -186,20 +286,20 @@ /*marshall DofIndexing data: */ - memcpy(marshalled_dataset,&numberofdofs,sizeof(numberofdofs));marshalled_dataset+=sizeof(numberofdofs); + memcpy(marshalled_dataset,&gsize,sizeof(gsize));marshalled_dataset+=sizeof(gsize); + memcpy(marshalled_dataset,&fsize,sizeof(fsize));marshalled_dataset+=sizeof(fsize); + memcpy(marshalled_dataset,&ssize,sizeof(ssize));marshalled_dataset+=sizeof(ssize); + memcpy(marshalled_dataset,&flagdoftype,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype); memcpy(marshalled_dataset,&clone,sizeof(clone));marshalled_dataset+=sizeof(clone); - memcpy(marshalled_dataset,f_set,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - memcpy(marshalled_dataset,s_set,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - memcpy(marshalled_dataset,doflist,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - - /*marshall doftype if there are some and add a flag*/ - if(this->doftype){ - flagdoftype=1; - memcpy(marshalled_dataset,&flagdoftype,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype); - memcpy(marshalled_dataset,doftype,numberofdofs*sizeof(int));marshalled_dataset+=numberofdofs*sizeof(int); - } - else{ - flagdoftype=0; - memcpy(marshalled_dataset,&flagdoftype,sizeof(flagdoftype));marshalled_dataset+=sizeof(flagdoftype); - } + + + if(this->gsize){ + memcpy(marshalled_dataset,f_set,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool); + memcpy(marshalled_dataset,s_set,gsize*sizeof(bool));marshalled_dataset+=gsize*sizeof(bool); + memcpy(marshalled_dataset,svalues,gsize*sizeof(double)); marshalled_dataset+=gsize*sizeof(double); + if(flagdoftype){ memcpy(marshalled_dataset,doftype,gsize*sizeof(int)); marshalled_dataset+=gsize*sizeof(int); } + memcpy(marshalled_dataset,gdoflist,gsize*sizeof(int)); marshalled_dataset+=gsize*sizeof(int); + } + if(this->fsize){ memcpy(marshalled_dataset,fdoflist,fsize*sizeof(int)); marshalled_dataset+=fsize*sizeof(int);} + if(this->ssize){ memcpy(marshalled_dataset,sdoflist,ssize*sizeof(int)); marshalled_dataset+=ssize*sizeof(int);} *pmarshalled_dataset=marshalled_dataset; @@ -210,19 +310,19 @@ int DofIndexing::MarshallSize(){ - int doftype_size=0; - - /*FInd size of doftype*/ - doftype_size+=sizeof(int); //Flag 0 or 1 - if (this->doftype){ - doftype_size+=numberofdofs*sizeof(int); - } - - return sizeof(numberofdofs)+ - sizeof(clone)+ - numberofdofs*sizeof(int)+ - numberofdofs*sizeof(int)+ - numberofdofs*sizeof(int)+ - doftype_size+ - sizeof(int); //sizeof(int) for enum type -} -/*}}}*/ + int size=0; + + size+=4*sizeof(int)+sizeof(bool); + if(this->gsize){ + size+= 2*this->gsize*sizeof(bool)+ + this->gsize*sizeof(double)+ + this->gsize*sizeof(int); + if(this->doftype)size+=this->gsize*sizeof(int); + } + if(this->fsize)size+=this->fsize*sizeof(int); + if(this->ssize)size+=this->ssize*sizeof(int); + + size+=sizeof(int); //sizeof(int) for enum type + + return size; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/DofIndexing.h =================================================================== --- /issm/trunk/src/c/objects/DofIndexing.h (revision 5771) +++ /issm/trunk/src/c/objects/DofIndexing.h (revision 5772) @@ -10,5 +10,8 @@ public: - int numberofdofs; //number of dofs for a node, variable + /*sizes: */ + int gsize; //number of dofs for a node + int fsize; //number of dofs solver for + int ssize; //number of constrained dofs /*partitioning: */ @@ -16,15 +19,22 @@ /*boundary conditions sets: */ - int* f_set; //set on which we solve - int* s_set; //set on which boundary conditions (dirichlet) are applied + bool* f_set; //is dof on f-set (on which we solve) + bool* s_set; //is dof on s-set (on which boundary conditions -dirichlet- are applied) + double* svalues; //list of constraint values. size g_size, for ease of use. + /*types of dofs: */ + int* doftype; //approximation type of the dofs (used only for coupling), size g_size + /*list of degrees of freedom: */ - int* doflist; //dof list on which we solve - int* doftype; //approximation type of the dofs (used only for coupling) + int* gdoflist; //dof list in g_set + int* fdoflist; //dof list in f_set + int* sdoflist; //dof list in s_set + /*DofIndexing constructors, destructors {{{1*/ DofIndexing(); - DofIndexing(int numberofdofs); - void Init(int numberofdofs,int* doftype); + DofIndexing(int g_size); + void Init(int g_size,int* doftype); + void InitSet(int setenum); DofIndexing(DofIndexing* properties); ~DofIndexing(); Index: /issm/trunk/src/c/objects/Elements/Element.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Element.h (revision 5771) +++ /issm/trunk/src/c/objects/Elements/Element.h (revision 5772) @@ -28,6 +28,6 @@ virtual void Configure(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0; virtual void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters)=0; - virtual void CreateKMatrix(Mat Kgg)=0; - virtual void CreatePVector(Vec pg)=0; + virtual void CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs)=0; + virtual void CreatePVector(Vec pg, Vec pf)=0; virtual void GetSolutionFromInputs(Vec solution)=0; virtual int GetNodeIndex(Node* node)=0; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 5772) @@ -367,5 +367,12 @@ /*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */ if (analysis_type==DiagnosticHorizAnalysisEnum){ - InputUpdateFromSolutionDiagnosticHoriz( solution); + int approximation; + inputs->GetParameterValue(&approximation,ApproximationEnum); + if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){ + InputUpdateFromSolutionDiagnosticStokes( solution); + } + else{ + InputUpdateFromSolutionDiagnosticHoriz( solution); + } } else if (analysis_type==DiagnosticHutterAnalysisEnum){ @@ -685,5 +692,5 @@ /*}}}*/ /*FUNCTION Penta::CreateKMatrix {{{1*/ -void Penta::CreateKMatrix(Mat Kgg){ +void Penta::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ int analysis_type; @@ -773,5 +780,5 @@ /*}}}*/ /*FUNCTION Penta::CreatePVector {{{1*/ -void Penta::CreatePVector(Vec pg){ +void Penta::CreatePVector(Vec pg, Vec pf){ int analysis_type; @@ -894,11 +901,7 @@ GetSolutionFromInputsDiagnosticStokes(solution); } - else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum ||approximation==HutterApproximationEnum){ + else{ GetSolutionFromInputsDiagnosticHoriz(solution); } - else if(approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum){ - return; //do nothing the elements around with update the solution - } - else ISSMERROR("approximation type : %i (%s) not implemented yet",approximation,EnumToString(approximation)); } else if(analysis_type==DiagnosticHutterAnalysisEnum){ @@ -1892,7 +1895,4 @@ this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum)); } - else if (*(iomodel->elements_type+index)==PattynStokesApproximationEnum){ - this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum)); - } else if (*(iomodel->elements_type+index)==NoneApproximationEnum){ this->inputs->AddInput(new IntInput(ApproximationEnum,NoneApproximationEnum)); @@ -1988,5 +1988,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); + tria->CreateKMatrix(Kgg,NULL,NULL); delete tria->matice; delete tria; @@ -2018,5 +2018,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); + tria->CreateKMatrix(Kgg,NULL,NULL); delete tria->matice; delete tria; @@ -2096,6 +2096,6 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICESp); - tria->GetDofList(&doflistm,MacAyealApproximationEnum); //Pattyn dof list - GetDofList(&doflistp,PattynApproximationEnum); //MacAyeal dof list + tria->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum); //Pattyn dof list + GetDofList(&doflistp,PattynApproximationEnum,GsetEnum); //MacAyeal dof list /*Retrieve all inputs we will be needing: */ @@ -2183,5 +2183,5 @@ if(IsOnWater())return; - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Spawn 3 beam elements: */ @@ -2321,5 +2321,5 @@ /*Call Tria function*/ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); + tria->CreateKMatrix(Kgg,NULL,NULL); delete tria->matice; delete tria; @@ -2343,5 +2343,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES); - tria->GetDofList(&doflist,MacAyealApproximationEnum); //Pattyn dof list + tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum); //Pattyn dof list /*Retrieve all inputs we will be needing: */ @@ -2399,5 +2399,5 @@ * grids: */ - tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg); + tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg,NULL,NULL); } @@ -2476,5 +2476,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,PattynApproximationEnum); + GetDofList(&doflist,PattynApproximationEnum,GsetEnum); /*Retrieve all inputs we will be needing: */ @@ -2596,5 +2596,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,StokesApproximationEnum); + GetDofList(&doflist,StokesApproximationEnum,GsetEnum); /*Retrieve all inputs we will be needing: */ @@ -2756,5 +2756,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /* Start looping on the number of gaussian points: */ @@ -2827,5 +2827,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); + tria->CreateKMatrix(Kgg,NULL,NULL); delete tria->matice; delete tria; @@ -2852,5 +2852,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreateKMatrix(Kgg); + tria->CreateKMatrix(Kgg,NULL,NULL); delete tria->matice; delete tria; return; @@ -2920,5 +2920,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); // /*recovre material parameters: */ @@ -3109,5 +3109,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); + tria->CreatePVector(pg,NULL); delete tria->matice; delete tria; @@ -3137,5 +3137,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); + tria->CreatePVector(pg,NULL); delete tria->matice; delete tria; @@ -3178,5 +3178,5 @@ /*recover doflist: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*recover parameters: */ @@ -3275,5 +3275,5 @@ *and use its CreatePVector functionality to return an elementary load vector: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); + tria->CreatePVector(pg,NULL); delete tria->matice; delete tria; return; @@ -3317,5 +3317,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,PattynApproximationEnum); + GetDofList(&doflist,PattynApproximationEnum,GsetEnum); /*Retrieve all inputs we will be needing: */ @@ -3429,5 +3429,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,StokesApproximationEnum); + GetDofList(&doflist,StokesApproximationEnum,GsetEnum); /*Retrieve all inputs we will be needing: */ @@ -3593,5 +3593,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Retrieve all inputs we will be needing: */ @@ -3656,5 +3656,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); + tria->CreatePVector(pg,NULL); delete tria->matice; delete tria; @@ -3680,5 +3680,5 @@ /*Spawn Tria element from the base of the Penta: */ tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria. - tria->CreatePVector(pg); + tria->CreatePVector(pg,NULL); delete tria->matice; delete tria; return; @@ -3744,5 +3744,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*recovre material parameters: */ @@ -3821,5 +3821,5 @@ /*}}}*/ /*FUNCTION Penta::GetDofList {{{1*/ -void Penta::GetDofList(int** pdoflist,int approximation_enum){ +void Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){ int i,j; @@ -3832,5 +3832,5 @@ /*First, figure out size of doflist: */ for(i=0;i<6;i++){ - numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum); + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); } @@ -3841,6 +3841,6 @@ count=0; for(i=0;i<6;i++){ - nodes[i]->GetDofList(doflist+count,approximation_enum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum); + nodes[i]->GetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); } @@ -3995,5 +3995,10 @@ /*If the element is a coupling, do nothing: every grid is also on an other elements * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/ - GetDofList(&doflist,approximation); + if(approximation==MacAyealPattynApproximationEnum){ + return; + } + else{ + GetDofList(&doflist,approximation,GsetEnum); + } /*Ok, we have vx and vy in values, fill in vx and vy arrays: */ @@ -4032,5 +4037,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -4070,5 +4075,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); Input* vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input); @@ -4106,5 +4111,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input); Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input); @@ -4152,5 +4157,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); Input* t_input=inputs->GetInput(TemperatureEnum); ISSMASSERT(t_input); @@ -4442,12 +4447,6 @@ InputUpdateFromSolutionDiagnosticPattyn(solution); } - else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){ - InputUpdateFromSolutionDiagnosticStokes(solution); - } else if (approximation==MacAyealPattynApproximationEnum){ InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution); - } - else if (approximation==PattynStokesApproximationEnum){ - InputUpdateFromSolutionDiagnosticPattynStokes(solution); } } @@ -4479,5 +4478,5 @@ /*Get dof list: */ - GetDofList(&doflist,MacAyealApproximationEnum); + GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -4572,6 +4571,6 @@ /*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */ - GetDofList(&doflistp,PattynApproximationEnum); - penta->GetDofList(&doflistm,MacAyealApproximationEnum); + GetDofList(&doflistp,PattynApproximationEnum,GsetEnum); + penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum); /*Get node data: */ @@ -4657,5 +4656,5 @@ /*Get dof list: */ - GetDofList(&doflist,PattynApproximationEnum); + GetDofList(&doflist,PattynApproximationEnum,GsetEnum); /*Get node data: */ @@ -4712,91 +4711,4 @@ } -/*}}}*/ -/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{1*/ -void Penta::InputUpdateFromSolutionDiagnosticPattynStokes(double* solution){ - - int i; - - const int numdofpervertexp=2; - const int numdofpervertexs=4; - const int numdofp=numdofpervertexp*NUMVERTICES; - const int numdofs=numdofpervertexs*NUMVERTICES; - int* doflistp=NULL; - int* doflists=NULL; - double pattyn_values[numdofp]; - double stokes_values[numdofs]; - double vx[NUMVERTICES]; - double vy[NUMVERTICES]; - double vz[NUMVERTICES]; - double vel[NUMVERTICES]; - int dummy; - double pressure[NUMVERTICES]; - double xyz_list[NUMVERTICES][3]; - - double *vz_ptr = NULL; - Penta *penta = NULL; - - /*OK, we have to add results of this element for pattyn - * and results from the penta at base for macayeal. Now recover results*/ - - /*Get dof listof this element (pattyn dofs and stokes dofs): */ - GetDofList(&doflistp,PattynApproximationEnum); - GetDofList(&doflists,StokesApproximationEnum); - - /*Get node data: */ - GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - - /*Use the dof list to index into the solution vector: */ - for(i=0;iGetInput(VzEnum); - if (vz_input){ - if (vz_input->Enum()!=PentaVertexInputEnum){ - ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum())); - } - vz_input->GetValuesPtr(&vz_ptr,&dummy); - for(i=0;iinputs->ChangeEnum(VxEnum,VxOldEnum); - this->inputs->ChangeEnum(VyEnum,VyOldEnum); - this->inputs->ChangeEnum(VzEnum,VzOldEnum); - this->inputs->ChangeEnum(PressureEnum,PressureOldEnum); - - /*Add vx and vy as inputs to the tria element: */ - this->inputs->AddInput(new PentaVertexInput(VxEnum,vx)); - this->inputs->AddInput(new PentaVertexInput(VyEnum,vy)); - this->inputs->AddInput(new PentaVertexInput(VzEnum,vz)); - this->inputs->AddInput(new PentaVertexInput(VelEnum,vel)); - this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure)); - - /*Free ressources:*/ - xfree((void**)&doflistp); - xfree((void**)&doflists); -} /*}}}*/ /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{1*/ @@ -4821,5 +4733,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get node data: */ @@ -4900,5 +4812,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get node data: */ @@ -4971,5 +4883,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5029,5 +4941,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5067,5 +4979,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5102,5 +5014,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5136,5 +5048,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5164,5 +5076,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector and extrude it */ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 5771) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 5772) @@ -75,6 +75,6 @@ void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); double RegularizeInversion(void); - void CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); void DeleteResults(void); int GetNodeIndex(Node* node); @@ -144,10 +144,10 @@ void CreatePVectorSlope( Vec pg); void CreatePVectorThermal( Vec pg); - void GetDofList(int** pdoflist,int approximation_enum=0); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetDofList1(int* doflist); - int GetElementType(void); - void GetParameterListOnVertices(double* pvalue,int enumtype); - void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); - void GetParameterValue(double* pvalue,Node* node,int enumtype); + int GetElementType(void); + void GetParameterListOnVertices(double* pvalue,int enumtype); + void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); + void GetParameterValue(double* pvalue,Node* node,int enumtype); void GetPhi(double* phi, double* epsilon, double viscosity); void GetSolutionFromInputsDiagnosticHoriz(Vec solutiong); Index: /issm/trunk/src/c/objects/Elements/Tria.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Elements/Tria.cpp (revision 5772) @@ -693,5 +693,5 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrix {{{1*/ -void Tria::CreateKMatrix(Mat Kgg){ +void Tria::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ /*retrive parameters: */ @@ -705,5 +705,5 @@ /*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */ if (analysis_type==DiagnosticHorizAnalysisEnum || analysis_type==AdjointHorizAnalysisEnum){ - CreateKMatrixDiagnosticMacAyeal( Kgg); + CreateKMatrixDiagnosticMacAyeal( Kgg,Kff,Kfs); } else if (analysis_type==DiagnosticHutterAnalysisEnum){ @@ -739,5 +739,5 @@ /*}}}*/ /*FUNCTION Tria::CreatePVector {{{1*/ -void Tria::CreatePVector(Vec pg){ +void Tria::CreatePVector(Vec pg, Vec pf){ int analysis_type; @@ -2440,5 +2440,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*retrieve some parameters: */ @@ -2578,5 +2578,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); this->parameters->FindParam(&dim,DimEnum); @@ -2680,5 +2680,5 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Retrieve all inputs we will be needed*/ @@ -2794,13 +2794,45 @@ /*}}}*/ /*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/ -void Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg){ +void Tria::CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs){ + + const int numdof = 2 *NUMVERTICES; + double *Ke_gg_viscous = NULL; + double *Ke_gg_friction = NULL; + ElementMatrix *Ke = NULL; + + /*compute stiffness matrices on the g-set, at the element level: */ + Ke_gg_viscous = CreateKMatrixDiagnosticMacAyealViscous(); + Ke_gg_friction = CreateKMatrixDiagnosticMacAyealFriction(); + + /*Initialize element matrix: */ + Ke=this->NewElementMatrix(MacAyealApproximationEnum); + + /*Add Ke_gg values to Ke element stifness matrix: */ + if(Ke_gg_viscous) Ke->AddValues(Ke_gg_viscous); + if(Ke_gg_friction)Ke->AddValues(Ke_gg_friction); + + /*Add Ke element stiffness matrix to global stiffness matrix: */ + Ke->AddToGlobal(Kgg,Kff,Kfs); + + /*Free ressources:*/ + delete Ke; + xfree((void**)&Ke_gg_viscous); + xfree((void**)&Ke_gg_friction); + +} + +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/ +double* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){ /* local declarations */ int i,j; + + /*output: */ + double* Ke_gg=NULL; /* node data: */ const int numdof=2*NUMVERTICES; double xyz_list[NUMVERTICES][3]; - int* doflist=NULL; /* gaussian points: */ @@ -2827,5 +2859,4 @@ /* local element matrices: */ - double Ke_gg[numdof][numdof]={0.0}; double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point. @@ -2834,4 +2865,5 @@ /*input parameters for structural analysis (diagnostic): */ double thickness; + /*retrieve some parameters: */ @@ -2839,9 +2871,11 @@ /*First, if we are on water, return empty matrix: */ - if(IsOnWater()) return; + if(IsOnWater()) return Ke_gg; + + /*allocate output: */ + Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double)); /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflist,MacAyealApproximationEnum); /*Retrieve all inputs we will be needing: */ @@ -2892,24 +2926,45 @@ /* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */ - for( i=0; iNewElementMatrix(MacAyealApproximationEnum); + + /*Add Ke_gg values to Ke element stifness matrix: */ + Ke->AddValues(Ke_gg_friction); + + /*Add Ke element stiffness matrix to global stiffness matrix: */ + Ke->AddToGlobal(Kgg,Kff,Kfs); + + /*Free ressources:*/ + delete Ke; + } + + /*Free ressources:*/ + xfree((void**)&Ke_gg_friction); +} + +/*}}}*/ +/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/ +double* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){ /* local declarations */ @@ -2917,9 +2972,10 @@ int analysis_type; + /*output: */ + double* Ke_gg=NULL; + /* node data: */ const int numdof = 2 *NUMVERTICES; double xyz_list[NUMVERTICES][3]; - int* doflistm=NULL; - int* doflistp=NULL; int numberofdofspernode=2; @@ -2934,5 +2990,4 @@ /* local element matrices: */ - double Ke_gg[numdof][numdof]={0.0}; double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag @@ -2952,4 +3007,5 @@ /*inputs: */ int drag_type; + /*retrive parameters: */ @@ -2965,11 +3021,12 @@ /* Get node coordinates and dof list: */ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); - GetDofList(&doflistm,MacAyealApproximationEnum); - GetDofList(&doflistp,PattynApproximationEnum); if (IsOnShelf()){ /*no friction, do nothing*/ - return; - } + return Ke_gg; + } + + /*allocte Ke_gg matrix: */ + Ke_gg=(double*)xcalloc(numdof*numdof,sizeof(double)); /*build friction object, used later on: */ @@ -3013,21 +3070,16 @@ &Ke_gg_gaussian[0][0],0); - for( i=0; iGetNumberOfDofs(approximation_enum); + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); } @@ -5251,6 +5307,6 @@ count=0; for(i=0;i<3;i++){ - nodes[i]->GetDofList(doflist+count,approximation_enum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum); + nodes[i]->GetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); } @@ -5258,4 +5314,97 @@ *pdoflist=doflist; +} +/*}}}*/ +/*FUNCTION Tria::GetInternalDofList {{{1*/ +int* Tria::GetInternalDofList(int approximation_enum,int setenum){ + + int i,j; + int numberofdofs=0; + int count=0; + int count2=0; + int ndof; + int ngdof; + + /*output: */ + int* doflist=NULL; + + /*First, figure out size of doflist: */ + for(i=0;i<3;i++){ + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + if(numberofdofs){ + /*Allocate: */ + doflist=(int*)xmalloc(numberofdofs*sizeof(int)); + + /*Populate: */ + count=0; + for(i=0;i<3;i++){ + nodes[i]->GetInternalDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + /*We now have something like: [0 1 0 2 0 1]. Offset by gsize, to get something like: + * [0 1 2 4 5 6]:*/ + count=0; + count2=0; + for(i=0;i<3;i++){ + ndof=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + for(j=count;jGetNumberOfDofs(approximation_enum,GsetEnum); + } + + } + else doflist=NULL; + + return doflist; + +} +/*}}}*/ +/*FUNCTION Tria::GetExternalDofList {{{1*/ +int* Tria::GetExternalDofList(int approximation_enum,int setenum){ + + int i,j; + int numberofdofs=0; + int count=0; + + /*output: */ + int* doflist=NULL; + + /*First, figure out size of doflist: */ + for(i=0;i<3;i++){ + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + if(numberofdofs){ + /*Allocate: */ + doflist=(int*)xmalloc(numberofdofs*sizeof(int)); + + /*Populate: */ + count=0; + for(i=0;i<3;i++){ + nodes[i]->GetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + } + else doflist=NULL; + + return doflist; +} +/*}}}*/ +/*FUNCTION Tria::GetNumberOfDofs {{{1*/ +int Tria::GetNumberOfDofs(int approximation_enum,int setenum){ + + int i; + + /*output: */ + int numberofdofs=0; + + for(i=0;i<3;i++){ + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + return numberofdofs; } /*}}}*/ @@ -5350,5 +5499,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get inputs*/ @@ -5393,5 +5542,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get inputs*/ @@ -5586,5 +5735,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5619,5 +5768,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5660,5 +5809,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5724,5 +5873,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5789,5 +5938,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Use the dof list to index into the solution vector: */ @@ -5830,4 +5979,97 @@ } /*}}}*/ +/*FUNCTION Tria::NewElementMatrix{{{1*/ +ElementMatrix* Tria::NewElementMatrix(int approximation){ + + /*parameters: */ + bool kff=false; + + /*output: */ + ElementMatrix* Ke=NULL; + int gsize; + int fsize; + int ssize; + int* gexternaldoflist=NULL; + int* finternaldoflist=NULL; + int* fexternaldoflist=NULL; + int* sinternaldoflist=NULL; + int* sexternaldoflist=NULL; + bool symmetric=true; + + /*retrieve some parameters: */ + this->parameters->FindParam(&kff,KffEnum); + + /*get number of dofs in sets g,f and s: */ + gsize=this->GetNumberOfDofs(approximation,GsetEnum); + if(kff){ + fsize=this->GetNumberOfDofs(approximation,FsetEnum); + ssize=this->GetNumberOfDofs(approximation,SsetEnum); + } + + /*get dof lists for f and s set: */ + if(!kff) gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum); + else{ + finternaldoflist=this->GetInternalDofList(approximation,FsetEnum); + fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum); + sinternaldoflist=this->GetInternalDofList(approximation,SsetEnum); + sexternaldoflist=this->GetExternalDofList(approximation,SsetEnum); + } + + /*Use symmetric constructor for ElementMatrix: */ + if(!kff) Ke=new ElementMatrix(gsize,symmetric,gexternaldoflist); + else Ke=new ElementMatrix(gsize,symmetric,finternaldoflist,fexternaldoflist,fsize,sinternaldoflist,sexternaldoflist,ssize); + + /*Free ressources and return:*/ + xfree((void**)&gexternaldoflist); + xfree((void**)&finternaldoflist); + xfree((void**)&fexternaldoflist); + xfree((void**)&sinternaldoflist); + xfree((void**)&sexternaldoflist); + + return Ke; +} +/*}}}*/ +/*FUNCTION Tria::NewElementVector{{{1*/ +ElementVector* Tria::NewElementVector(int approximation){ + + /*parameters: */ + bool kff=false; + + /*output: */ + ElementVector* pe=NULL; + int gsize; + int fsize; + int* gexternaldoflist=NULL; + int* finternaldoflist=NULL; + int* fexternaldoflist=NULL; + + /*retrieve some parameters: */ + this->parameters->FindParam(&kff,KffEnum); + + /*get number of dofs in sets g,f and s: */ + gsize=this->GetNumberOfDofs(approximation,GsetEnum); + if(kff)fsize=this->GetNumberOfDofs(approximation,FsetEnum); + + /*get dof lists for f and s set: */ + if(!kff){ + gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum); + } + else{ + finternaldoflist=this->GetInternalDofList(approximation,FsetEnum); + fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum); + } + + /*Use symmetric constructor for ElementVector: */ + if(!kff)pe=new ElementVector(gsize,gexternaldoflist); + else pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize); + + /*Free ressources and return:*/ + xfree((void**)&gexternaldoflist); + xfree((void**)&finternaldoflist); + xfree((void**)&fexternaldoflist); + + return pe; +} +/*}}}*/ /*FUNCTION Tria::SetClone {{{1*/ void Tria::SetClone(int* minranks){ Index: /issm/trunk/src/c/objects/Elements/Tria.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Tria.h (revision 5771) +++ /issm/trunk/src/c/objects/Elements/Tria.h (revision 5772) @@ -17,4 +17,6 @@ class Matice; class Matpar; +class ElementMatrix; +class ElementVector; #include "../../include/include.h" @@ -72,6 +74,6 @@ void SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters); double RegularizeInversion(void); - void CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); int GetNodeIndex(Node* node); bool IsOnBed(); @@ -122,7 +124,9 @@ void CreateKMatrixBalancedthickness_CG(Mat Kgg); void CreateKMatrixBalancedvelocities(Mat Kgg); - void CreateKMatrixDiagnosticMacAyeal(Mat Kgg); + void CreateKMatrixDiagnosticMacAyeal(Mat Kgg,Mat Kff, Mat Kfs); + double* CreateKMatrixDiagnosticMacAyealViscous(void); + double* CreateKMatrixDiagnosticMacAyealFriction(void); + void CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg,Mat Kff, Mat Kfs); void CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg); - void CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg); void CreateKMatrixDiagnosticPattynFriction(Mat Kgg); void CreateKMatrixDiagnosticHutter(Mat Kgg); @@ -147,14 +151,17 @@ void CreatePVectorThermalSheet( Vec pg); void CreatePVectorThermalShelf( Vec pg); - double GetArea(void); - int GetElementType(void); - void GetDofList(int** pdoflist,int approximation_enum=0); + double GetArea(void); + int GetElementType(void); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetDofList1(int* doflist); - void GetParameterListOnVertices(double* pvalue,int enumtype); - void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); - void GetParameterValue(double* pvalue,Node* node,int enumtype); + int* GetInternalDofList(int approximation_enum,int setenum); + int* GetExternalDofList(int approximation_enum,int setenum); + int GetNumberOfDofs(int approximation_enum,int setenum); + void GetParameterListOnVertices(double* pvalue,int enumtype); + void GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue); + void GetParameterValue(double* pvalue,Node* node,int enumtype); void GetSolutionFromInputsDiagnosticHoriz(Vec solution); void GetSolutionFromInputsDiagnosticHutter(Vec solution); - void GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input); + void GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input); void GradjDragStokes(Vec gradient); void InputUpdateFromSolutionAdjointBalancedthickness( double* solution); @@ -164,4 +171,6 @@ void InputUpdateFromSolutionOneDof(double* solution,int enum_type); bool IsInput(int name); + ElementMatrix* NewElementMatrix(int approximation); + ElementVector* NewElementVector(int approximation); void SetClone(int* minranks); void SurfaceNormal(double* surface_normal, double xyz_list[3][3]); Index: /issm/trunk/src/c/objects/FemModel.cpp =================================================================== --- /issm/trunk/src/c/objects/FemModel.cpp (revision 5771) +++ /issm/trunk/src/c/objects/FemModel.cpp (revision 5772) @@ -38,5 +38,4 @@ analysis_type_list=(int*)xmalloc(nummodels*sizeof(int)); m_nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*)); - m_yg=(Vec*)xmalloc(nummodels*sizeof(Vec)); m_ys=(Vec*)xmalloc(nummodels*sizeof(Vec)); @@ -44,5 +43,4 @@ for(i=0;iSetCurrentConfiguration(analysis_type); - _printf_(" create degrees of freedom\n"); + _printf_(" create vertex degrees of freedom\n"); VerticesDofx(&partition,&tpartition,vertices,parameters); + + _printf_(" resolve node constraints\n"); + SpcNodesx(nodes,constraints,analysis_type); + + _printf_(" create nodal degrees of freedom\n"); NodesDofx(nodes,parameters,analysis_type); - - _printf_(" create single point constraints\n"); - SpcNodesx( &m_yg[i], nodes,constraints,analysis_type); + + _printf_(" create nodal constraints vector\n"); + CreateNodalConstraintsx(&m_ys[i],nodes,analysis_type); _printf_(" create node sets\n"); BuildNodeSetsx(&m_nodesets[i], nodes,analysis_type); - - _printf_(" reducing single point constraints vector\n"); - Reducevectorgtosx(&m_ys[i], m_yg[i],m_nodesets[i]); _printf_(" configuring element and loads\n"); @@ -97,6 +97,4 @@ NodeSets* temp_nodesets=m_nodesets[i]; delete temp_nodesets; - Vec temp_yg=m_yg[i]; - VecFree(&temp_yg); Vec temp_ys=m_ys[i]; VecFree(&temp_ys); @@ -105,5 +103,4 @@ /*Delete dynamically allocated arrays: */ xfree((void**)&m_nodesets); - xfree((void**)&m_yg); xfree((void**)&m_ys); @@ -147,5 +144,4 @@ /*activate matrices/vectors: */ nodesets=m_nodesets[analysis_counter]; - yg=m_yg[analysis_counter]; ys=m_ys[analysis_counter]; Index: /issm/trunk/src/c/objects/FemModel.h =================================================================== --- /issm/trunk/src/c/objects/FemModel.h (revision 5771) +++ /issm/trunk/src/c/objects/FemModel.h (revision 5772) @@ -41,5 +41,4 @@ //multiple sets of matrices/vectors for each analysis_type. m stands for multiple NodeSets** m_nodesets; //boundary conditions dof sets - Vec* m_yg; //boundary conditions in global g-set Vec* m_ys; //boundary conditions, in reduced s-set Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 5771) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 5772) @@ -188,4 +188,5 @@ IoModelFetchData(&this->stokesreconditioning,iomodel_handle,"stokesreconditioning"); IoModelFetchData(&this->waitonlock,iomodel_handle,"waitonlock"); + IoModelFetchData(&this->kff,iomodel_handle,"kff"); /*!Get thermal parameters: */ @@ -324,4 +325,5 @@ this->connectivity=0; this->lowmem=0; + this->kff=0; this->optscal=NULL; this->yts=0; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 5771) +++ /issm/trunk/src/c/objects/IoModel.h (revision 5772) @@ -150,4 +150,5 @@ double yts; double waitonlock; + int kff; /*thermal parameters: */ Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Icefront.cpp (revision 5772) @@ -307,5 +307,5 @@ /*}}}*/ /*FUNCTION Icefront::CreateKMatrix {{{1*/ -void Icefront::CreateKMatrix(Mat Kgg){ +void Icefront::CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs){ /*No stiffness loads applied, do nothing: */ @@ -315,5 +315,5 @@ /*}}}*/ /*FUNCTION Icefront::CreatePVector {{{1*/ -void Icefront::CreatePVector(Vec pg){ +void Icefront::CreatePVector(Vec pg, Vec pf){ /*Checks in debugging mode*/ @@ -331,5 +331,5 @@ inputs->GetParameterValue(&type,TypeEnum); if (type==MacAyeal2dIceFrontEnum){ - CreatePVectorDiagnosticMacAyeal2d( pg); + CreatePVectorDiagnosticMacAyeal2d( pg,pf); } else if (type==MacAyeal3dIceFrontEnum){ @@ -353,10 +353,10 @@ /*}}}*/ /*FUNCTION Icefront::PenaltyCreateKMatrix {{{1*/ -void Icefront::PenaltyCreateKMatrix(Mat Kgg,double kmax){ +void Icefront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs, double kmax){ /*do nothing: */ } /*}}}*/ /*FUNCTION Icefront::PenaltyCreatePVector{{{1*/ -void Icefront::PenaltyCreatePVector(Vec pg,double kmax){ +void Icefront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ /*do nothing: */ } @@ -423,9 +423,12 @@ /*Icefront numerics: */ /*FUNCTION Icefront::CreatePVectorDiagnosticMacAyeal2d{{{1*/ -void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg){ +void Icefront::CreatePVectorDiagnosticMacAyeal2d( Vec pg,Vec pf){ /*Constants*/ const int numnodes= NUMVERTICESSEG; const int numdofs = numnodes *NDOF2; + + /*output: */ + ElementVector* pe=NULL; /*Intermediary*/ @@ -434,5 +437,4 @@ double thickness,bed,pressure,ice_pressure,rho_water,rho_ice,gravity; double water_pressure,air_pressure,surface_under_water,base_under_water; - int *doflist = NULL; double xyz_list[numnodes][3]; double pe_g[numdofs] = {0.0}; @@ -444,5 +446,4 @@ /* Get dof list and node coordinates: */ - GetDofList(&doflist,MacAyealApproximationEnum); GetVerticesCoordinates(&xyz_list[0][0],nodes,numnodes); @@ -493,10 +494,19 @@ } } - - VecSetValues(pg,numdofs,doflist,pe_g,ADD_VALUES); + + /*Initialize element matrix: */ + pe=this->NewElementVector(MacAyealApproximationEnum); + + /*Add pe_g values to pe element stifness load: */ + pe->AddValues(&pe_g[0]); + + /*Add pe element load vector to global load vector: */ + pe->AddToGlobal(pg,pf); + + /*Free ressources:*/ + delete pe; /*Free ressources:*/ delete gauss; - xfree((void**)&doflist); } /*}}}*/ @@ -521,5 +531,5 @@ icefront->element=tria; icefront->inputs->AddInput(new IntInput(TypeEnum,MacAyeal2dIceFrontEnum)); - icefront->CreatePVectorDiagnosticMacAyeal2d(pg); + icefront->CreatePVectorDiagnosticMacAyeal2d(pg,NULL); /*clean-up*/ @@ -574,5 +584,5 @@ /* Get dof list and node coordinates: */ - GetDofList(&doflist,PattynApproximationEnum); + GetDofList(&doflist,PattynApproximationEnum,GsetEnum); GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); @@ -692,5 +702,5 @@ /* Get dof list and node coordinates: */ - GetDofList(&doflist,StokesApproximationEnum); + GetDofList(&doflist,StokesApproximationEnum,GsetEnum); GetVerticesCoordinates(&xyz_list[0][0], nodes, numnodes); @@ -766,5 +776,5 @@ /*}}}*/ /*FUNCTION Icefront::GetDofList {{{1*/ -void Icefront::GetDofList(int** pdoflist,int approximation_enum){ +void Icefront::GetDofList(int** pdoflist,int approximation_enum,int setenum){ int i,j; @@ -792,5 +802,5 @@ /*Figure out size of doflist: */ for(i=0;iGetNumberOfDofs(approximation_enum); + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); } @@ -801,10 +811,140 @@ count=0; for(i=0;iGetDofList(doflist+count,approximation_enum); - count+=nodes[i]->GetNumberOfDofs(approximation_enum); + nodes[i]->GetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); } /*Assign output pointers:*/ *pdoflist=doflist; +} +/*}}}*/ +/*FUNCTION Icefront::GetInternalDofList {{{1*/ +int* Icefront::GetInternalDofList(int approximation_enum,int setenum){ + + int i,j; + int numberofdofs=0; + int count=0; + int count2=0; + int type; + int numberofnodes=2; + int ndof; + + /*output: */ + int* doflist=NULL; + + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + /*Some checks for debugging*/ + ISSMASSERT(nodes); + + /*How many nodes? :*/ + if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum) + numberofnodes=2; + else + numberofnodes=4; + + /*Figure out size of doflist: */ + for(i=0;iGetNumberOfDofs(approximation_enum,setenum); + } + + if(numberofdofs){ + /*Allocate: */ + doflist=(int*)xmalloc(numberofdofs*sizeof(int)); + + /*Populate: */ + count=0; + for(i=0;iGetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + + /*We now have something like: [0 1 0 1]. Offset by gsize, to get something like: + * [0 1 2 3]: */ + count=0; + count2=0; + for(i=0;iGetNumberOfDofs(approximation_enum,setenum); + for(j=count;jGetNumberOfDofs(approximation_enum,GsetEnum); + } + } + + return doflist; +} +/*}}}*/ +/*FUNCTION Icefront::GetExternalDofList{{{1*/ +int* Icefront::GetExternalDofList(int approximation_enum,int setenum){ + + int i,j; + int numberofdofs=0; + int count=0; + int type; + int numberofnodes=2; + + /*output: */ + int* doflist=NULL; + + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + /*Some checks for debugging*/ + ISSMASSERT(nodes); + + /*How many nodes? :*/ + if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum) + numberofnodes=2; + else + numberofnodes=4; + + /*Figure out size of doflist: */ + for(i=0;iGetNumberOfDofs(approximation_enum,setenum); + } + + if(numberofdofs){ + /*Allocate: */ + doflist=(int*)xmalloc(numberofdofs*sizeof(int)); + + /*Populate: */ + count=0; + for(i=0;iGetDofList(doflist+count,approximation_enum,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum); + } + } + + return doflist; +} +/*}}}*/ +/*FUNCTION Icefront::GetNumberOfDofs {{{1*/ +int Icefront::GetNumberOfDofs(int approximation_enum,int setenum){ + + + int i; + int numberofdofs=0; + int type; + int numberofnodes=2; + + /*recover type: */ + inputs->GetParameterValue(&type,TypeEnum); + + /*Some checks for debugging*/ + ISSMASSERT(nodes); + + /*How many nodes? :*/ + if(type==MacAyeal2dIceFrontEnum || type==MacAyeal3dIceFrontEnum) + numberofnodes=2; + else + numberofnodes=4; + + /*Figure out size of doflist: */ + for(i=0;iGetNumberOfDofs(approximation_enum,setenum); + } + + return numberofdofs; } /*}}}*/ @@ -1309,2 +1449,44 @@ } /*}}}*/ +/*FUNCTION Icefront::NewElementVector{{{1*/ +ElementVector* Icefront::NewElementVector(int approximation){ + + /*parameters: */ + bool kff=false; + + /*output: */ + ElementVector* pe=NULL; + int gsize; + int fsize; + int* gexternaldoflist=NULL; + int* finternaldoflist=NULL; + int* fexternaldoflist=NULL; + + /*retrieve some parameters: */ + this->parameters->FindParam(&kff,KffEnum); + + /*get number of dofs in sets g and f: */ + gsize=this->GetNumberOfDofs(approximation,GsetEnum); + if(kff)fsize=this->GetNumberOfDofs(approximation,FsetEnum); + + /*get dof lists for f and s set: */ + if(!kff){ + gexternaldoflist=this->GetExternalDofList(approximation,GsetEnum); + } + else{ + finternaldoflist=this->GetInternalDofList(approximation,FsetEnum); + fexternaldoflist=this->GetExternalDofList(approximation,FsetEnum); + } + + /*Use symmetric constructor for ElementVector: */ + if(!kff)pe=new ElementVector(gsize,gexternaldoflist); + else pe=new ElementVector(gsize,finternaldoflist,fexternaldoflist,fsize); + + /*Free ressources and return:*/ + xfree((void**)&gexternaldoflist); + xfree((void**)&finternaldoflist); + xfree((void**)&fexternaldoflist); + + return pe; +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Icefront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Icefront.h (revision 5772) @@ -16,4 +16,5 @@ class Element; class IoModel; +class ElementVector; /*}}}*/ @@ -69,16 +70,19 @@ 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 CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); - void PenaltyCreateKMatrix(Mat Kgg,double kmax); - void PenaltyCreatePVector(Vec pg,double kmax); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); + void PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax); + void PenaltyCreatePVector(Vec pg,Vec pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ /*Load management: {{{1*/ - void CreatePVectorDiagnosticMacAyeal2d(Vec pg); + void CreatePVectorDiagnosticMacAyeal2d(Vec pg,Vec pf); void CreatePVectorDiagnosticMacAyeal3d(Vec pg); void CreatePVectorDiagnosticPattyn(Vec pg); void CreatePVectorDiagnosticStokes(Vec pg); - void GetDofList(int** pdoflist,int approximation_enum=0); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); + int* GetInternalDofList(int approximation_enum,int setenum); + int* GetExternalDofList(int approximation_enum,int setenum); + int GetNumberOfDofs(int approximation_enum,int setenum); void QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); @@ -86,4 +90,5 @@ double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list); void GetNormal(double* normal,double xyz_list[2][3]); + ElementVector* NewElementVector(int approximation); /*}}}*/ }; Index: /issm/trunk/src/c/objects/Loads/Load.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Load.h (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Load.h (revision 5772) @@ -26,8 +26,8 @@ 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,Vertices* vertices,Materials* materials,Parameters* parameters)=0; - virtual void CreateKMatrix(Mat Kgg)=0; - virtual void CreatePVector(Vec pg)=0; - virtual void PenaltyCreateKMatrix(Mat Kgg,double kmax)=0; - virtual void PenaltyCreatePVector(Vec pg,double kmax)=0; + virtual void CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs)=0; + virtual void CreatePVector(Vec pg, Vec pf)=0; + virtual void PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs, double kmax)=0; + virtual void PenaltyCreatePVector(Vec pg, Vec pf, double kmax)=0; virtual bool InAnalysis(int analysis_type)=0; /*}}}*/ Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp (revision 5772) @@ -328,5 +328,5 @@ /*}}}*/ /*FUNCTION Numericalflux::CreateKMatrix {{{1*/ -void Numericalflux::CreateKMatrix(Mat Kgg){ +void Numericalflux::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){ int type; @@ -349,5 +349,5 @@ /*}}}*/ /*FUNCTION Numericalflux::CreatePVector {{{1*/ -void Numericalflux::CreatePVector(Vec pg){ +void Numericalflux::CreatePVector(Vec pg,Vec pf){ int type; @@ -374,5 +374,5 @@ /*}}}*/ /*FUNCTION Numericalflux::PenaltyCreateKMatrix {{{1*/ -void Numericalflux::PenaltyCreateKMatrix(Mat Kgg,double kmax){ +void Numericalflux::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ /*No stiffness loads applied, do nothing: */ @@ -382,5 +382,5 @@ /*}}}*/ /*FUNCTION Numericalflux::PenaltyCreatePVector{{{1*/ -void Numericalflux::PenaltyCreatePVector(Vec pg,double kmax){ +void Numericalflux::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ /*No penalty loads applied, do nothing: */ @@ -417,5 +417,5 @@ GetVerticesCoordinates(&xyz_list[0][0], nodes,NUMVERTICES_INTERNAL); - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Retrieve all inputs and parameters we will be needing: */ @@ -523,5 +523,5 @@ } - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /* Start looping on the number of gaussian points: */ @@ -617,5 +617,5 @@ } - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /* Start looping on the number of gaussian points: */ @@ -646,5 +646,5 @@ /*}}}*/ /*FUNCTION Numericalflux::GetDofList {{{1*/ -void Numericalflux::GetDofList(int** pdoflist){ +void Numericalflux::GetDofList(int** pdoflist,int approximation,int setenum){ int i,j; @@ -671,5 +671,5 @@ /*Figure out size of doflist: */ for(i=0;iGetNumberOfDofs(); + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum); } @@ -680,6 +680,6 @@ count=0; for(i=0;iGetDofList(doflist+count); - count+=nodes[i]->GetNumberOfDofs(); + nodes[i]->GetDofList(doflist+count,approximation,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation,setenum); } Index: /issm/trunk/src/c/objects/Loads/Numericalflux.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Numericalflux.h (revision 5772) @@ -65,12 +65,12 @@ 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 CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); - void PenaltyCreateKMatrix(Mat Kgg,double kmax); - void PenaltyCreatePVector(Vec pg,double kmax); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); + void PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax); + void PenaltyCreatePVector(Vec pg,Vec pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ /*Numericalflux management:{{{1*/ - void GetDofList(int** pdoflist); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); void GetNormal(double* normal,double xyz_list[4][3]); void CreateKMatrixInternal(Mat Kgg); Index: /issm/trunk/src/c/objects/Loads/Pengrid.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Pengrid.cpp (revision 5772) @@ -281,5 +281,5 @@ /*}}}1*/ /*FUNCTION Pengrid::CreateKMatrix {{{1*/ -void Pengrid::CreateKMatrix(Mat Kgg){ +void Pengrid::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){ /*No loads applied, do nothing: */ @@ -289,5 +289,5 @@ /*}}}1*/ /*FUNCTION Pengrid::CreatePVector {{{1*/ -void Pengrid::CreatePVector(Vec pg){ +void Pengrid::CreatePVector(Vec pg,Vec pf){ /*No loads applied, do nothing: */ @@ -297,5 +297,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreateMatrix {{{1*/ -void Pengrid::PenaltyCreateKMatrix(Mat Kgg,double kmax){ +void Pengrid::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ int analysis_type; @@ -320,5 +320,5 @@ /*}}}1*/ /*FUNCTION Pengrid::PenaltyCreatePVector {{{1*/ -void Pengrid::PenaltyCreatePVector(Vec pg,double kmax){ +void Pengrid::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ int analysis_type; @@ -419,5 +419,5 @@ /*Pengrid management:*/ /*FUNCTION Pengrid::GetDofList {{{1*/ -void Pengrid::GetDofList(int** pdoflist){ +void Pengrid::GetDofList(int** pdoflist,int approximation,int setenum){ int i,j; @@ -431,5 +431,5 @@ /*First, figure out size of doflist: */ - numberofdofs=node->GetNumberOfDofs(); + numberofdofs=node->GetNumberOfDofs(approximation,setenum); /*Allocate: */ @@ -437,5 +437,5 @@ /*Populate: */ - node->GetDofList(doflist); + node->GetDofList(doflist,approximation,setenum); /*Assign output pointers:*/ @@ -567,5 +567,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); //recover slope: */ @@ -620,5 +620,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); //Compute pressure melting point @@ -656,5 +656,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); Ke[0][0]=kmax*pow((double)10,penalty_offset); @@ -693,5 +693,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); //First recover pressure and temperature values, using the element: */ @@ -756,5 +756,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); //First recover pressure and penalty_offset Index: /issm/trunk/src/c/objects/Loads/Pengrid.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Pengrid.h (revision 5772) @@ -71,13 +71,13 @@ 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 CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); - void PenaltyCreateKMatrix(Mat Kgg,double kmax); - void PenaltyCreatePVector(Vec pg,double kmax); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); + void PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax); + void PenaltyCreatePVector(Vec pg,Vec pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ /*Pengrid management {{{1*/ void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax); - void GetDofList(int** pdoflist); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); void PenaltyCreateKMatrixThermal(Mat Kgg,double kmax); void PenaltyCreateKMatrixMelting(Mat Kgg,double kmax); Index: /issm/trunk/src/c/objects/Loads/Penpair.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Penpair.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Penpair.cpp (revision 5772) @@ -184,5 +184,5 @@ /*}}}1*/ /*FUNCTION Penpair::CreateKMatrix {{{1*/ -void Penpair::CreateKMatrix(Mat Kgg){ +void Penpair::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){ /*If you code this piece, don't forget that a penalty will be inactive if it is dealing with clone nodes*/ /*No loads applied, do nothing: */ @@ -192,5 +192,5 @@ /*}}}1*/ /*FUNCTION Penpair::CreatePVector {{{1*/ -void Penpair::CreatePVector(Vec pg){ +void Penpair::CreatePVector(Vec pg,Vec pf){ /*No loads applied, do nothing: */ @@ -200,5 +200,5 @@ /*}}}1*/ /*FUNCTION Penpair::PenaltyCreateKMatrix {{{1*/ -void Penpair::PenaltyCreateKMatrix(Mat Kgg,double kmax){ +void Penpair::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ int analysis_type; @@ -251,5 +251,5 @@ /*}}}1*/ /*FUNCTION Penpair::PenaltyCreatePVector {{{1*/ -void Penpair::PenaltyCreatePVector(Vec pg,double kmax){ +void Penpair::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ /*No loads applied, do nothing: */ return; @@ -297,5 +297,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*recover pointers: */ @@ -340,5 +340,5 @@ /*Get dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*recover pointers: */ @@ -378,5 +378,5 @@ /*}}}1*/ /*FUNCTION Penpair::GetDofList {{{1*/ -void Penpair::GetDofList(int** pdoflist){ +void Penpair::GetDofList(int** pdoflist,int approximation,int setenum){ int i,j; @@ -398,5 +398,5 @@ /*First, figure out size of doflist: */ for(i=0;i<2;i++){ - numberofdofs+=nodes[i]->GetNumberOfDofs(); + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum); } @@ -407,6 +407,6 @@ count=0; for(i=0;i<2;i++){ - nodes[i]->GetDofList(doflist+count); - count+=nodes[i]->GetNumberOfDofs(); + nodes[i]->GetDofList(doflist+count,approximation,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation,setenum); } Index: /issm/trunk/src/c/objects/Loads/Penpair.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Penpair.h (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Penpair.h (revision 5772) @@ -57,8 +57,8 @@ 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 CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); - void PenaltyCreateKMatrix(Mat Kgg,double kmax); - void PenaltyCreatePVector(Vec pg,double kmax); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); + void PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax); + void PenaltyCreatePVector(Vec pg,Vec pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ @@ -66,5 +66,5 @@ void PenaltyCreateKMatrixDiagnosticHoriz(Mat Kgg,double kmax); void PenaltyCreateKMatrixDiagnosticStokes(Mat Kgg,double kmax); - void GetDofList(int** pdoflist); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); /*}}}*/ }; Index: /issm/trunk/src/c/objects/Loads/Riftfront.cpp =================================================================== --- /issm/trunk/src/c/objects/Loads/Riftfront.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Riftfront.cpp (revision 5772) @@ -376,15 +376,15 @@ /*}}}*/ /*FUNCTION Riftfront::CreateKMatrix {{{1*/ -void Riftfront::CreateKMatrix(Mat Kgg){ +void Riftfront::CreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs){ /*do nothing: */ } /*}}}1*/ /*FUNCTION Riftfront::CreatePVector {{{1*/ -void Riftfront::CreatePVector(Vec pg){ +void Riftfront::CreatePVector(Vec pg,Vec pf){ /*do nothing: */ } /*}}}1*/ /*FUNCTION Riftfront::PenaltyCreateKMatrix {{{1*/ -void Riftfront::PenaltyCreateKMatrix(Mat Kgg,double kmax){ +void Riftfront::PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat Kfs,double kmax){ int i; @@ -419,5 +419,5 @@ /* Get node coordinates and dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /* Set Ke_gg to 0: */ @@ -499,5 +499,5 @@ /*}}}1*/ /*FUNCTION Riftfront::PenaltyCreatePVector {{{1*/ -void Riftfront::PenaltyCreatePVector(Vec pg,double kmax){ +void Riftfront::PenaltyCreatePVector(Vec pg,Vec pf,double kmax){ int i ,j; @@ -544,5 +544,5 @@ /* Get node coordinates and dof list: */ - GetDofList(&doflist); + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); /*Get some inputs: */ @@ -731,5 +731,5 @@ /*}}}1*/ /*FUNCTION Riftfront::GetDofList {{{1*/ -void Riftfront::GetDofList(int** pdoflist){ +void Riftfront::GetDofList(int** pdoflist,int approximation,int setenum){ int i,j; @@ -751,5 +751,5 @@ /*Figure out size of doflist: */ for(i=0;iGetNumberOfDofs(); + numberofdofs+=nodes[i]->GetNumberOfDofs(approximation,setenum); } @@ -760,6 +760,6 @@ count=0; for(i=0;iGetDofList(doflist+count); - count+=nodes[i]->GetNumberOfDofs(); + nodes[i]->GetDofList(doflist+count,approximation,setenum); + count+=nodes[i]->GetNumberOfDofs(approximation,setenum); } Index: /issm/trunk/src/c/objects/Loads/Riftfront.h =================================================================== --- /issm/trunk/src/c/objects/Loads/Riftfront.h (revision 5771) +++ /issm/trunk/src/c/objects/Loads/Riftfront.h (revision 5772) @@ -74,12 +74,12 @@ 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 CreateKMatrix(Mat Kgg); - void CreatePVector(Vec pg); - void PenaltyCreateKMatrix(Mat Kgg,double kmax); - void PenaltyCreatePVector(Vec pg,double kmax); + void CreateKMatrix(Mat Kgg, Mat Kff, Mat Kfs); + void CreatePVector(Vec pg, Vec pf); + void PenaltyCreateKMatrix(Mat Kgg,Mat Kff, Mat kfs, double kmax); + void PenaltyCreatePVector(Vec pg,Vec pf, double kmax); bool InAnalysis(int analysis_type); /*}}}*/ /*Riftfront specific routines: {{{1*/ - void GetDofList(int** pdoflist); + void GetDofList(int** pdoflist,int approximation_enum,int setenum); bool PreStable(); void SetPreStable(); Index: /issm/trunk/src/c/objects/Node.cpp =================================================================== --- /issm/trunk/src/c/objects/Node.cpp (revision 5771) +++ /issm/trunk/src/c/objects/Node.cpp (revision 5772) @@ -63,5 +63,5 @@ /*Intermediary*/ int k; - int numdofs; + int gsize; /*id: */ @@ -72,5 +72,5 @@ /*indexing:*/ DistributeNumDofs(&this->indexing,analysis_type,iomodel->vertices_type+io_index); //number of dofs per node - numdofs=this->indexing.numberofdofs; + gsize=this->indexing.gsize; /*Hooks*/ @@ -95,9 +95,9 @@ if (!iomodel->vertices_type) ISSMERROR("iomodel->vertices_type is NULL"); if (iomodel->vertices_type[io_index]==MacAyealApproximationEnum && !iomodel->gridonbed[io_index]){ - for(k=1;k<=numdofs;k++) this->FreezeDof(k); + for(k=1;k<=gsize;k++) this->FreezeDof(k); } if (iomodel->vertices_type[io_index]==MacAyealPattynApproximationEnum && iomodel->gridonmacayeal[io_index]){ if(!iomodel->gridonbed[io_index]){ - for(k=1;k<=numdofs;k++) this->FreezeDof(k); + for(k=1;k<=gsize;k++) this->FreezeDof(k); } } @@ -106,5 +106,5 @@ if (!iomodel->gridonhutter) ISSMERROR("iomodel->gridonhutter is NULL"); if (iomodel->gridonhutter[io_index]){ - for(k=1;k<=numdofs;k++){ + for(k=1;k<=gsize;k++){ this->FreezeDof(k); } @@ -117,5 +117,5 @@ if (!iomodel->gridonhutter) ISSMERROR("iomodel->gridonhutter is NULL"); if (!iomodel->gridonhutter[io_index]){ - for(k=1;k<=numdofs;k++){ + for(k=1;k<=gsize;k++){ this->FreezeDof(k); } @@ -136,5 +136,5 @@ if (!iomodel->gridonbed) ISSMERROR("iomodel->gridonbed is NULL"); if (!iomodel->gridonbed[io_index]){ - for(k=1;k<=numdofs;k++){ + for(k=1;k<=gsize;k++){ this->FreezeDof(k); } @@ -199,5 +199,5 @@ int enum_type=0; char* marshalled_inputs=NULL; - int marshalled_inputs_size; + int marshalled_inputssize; /*recover marshalled_dataset: */ @@ -220,8 +220,8 @@ /*Marshall inputs: */ - marshalled_inputs_size=inputs->MarshallSize(); + marshalled_inputssize=inputs->MarshallSize(); marshalled_inputs=inputs->Marshall(); - memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char)); - marshalled_dataset+=marshalled_inputs_size; + memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputssize*sizeof(char)); + marshalled_dataset+=marshalled_inputssize; /*Free ressources:*/ @@ -294,5 +294,5 @@ /*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective - * datasets, using internal ids and offsets hidden in hooks: */ + * datasets, using internal ids and off_sets hidden in hooks: */ hvertex->configure(verticesin); @@ -303,7 +303,16 @@ } /*FUNCTION Node::GetDof {{{1*/ -int Node::GetDof(int dofindex){ - - return indexing.doflist[dofindex]; +int Node::GetDof(int dofindex,int setenum){ + + if(setenum==GsetEnum){ + return indexing.gdoflist[dofindex]; + } + else if(setenum==FsetEnum){ + return indexing.fdoflist[dofindex]; + } + else if(setenum==SsetEnum){ + return indexing.sdoflist[dofindex]; + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); } @@ -320,21 +329,150 @@ /*}}}*/ /*FUNCTION Node::GetDofList{{{1*/ -void Node::GetDofList(int* outdoflist,int approximation_enum){ +void Node::GetDofList(int* outdoflist,int approximation_enum,int setenum){ int i; int count=0; - - if(approximation_enum && this->indexing.doftype){ //if an approximation is specified and this node has several approsimations - for(i=0;iindexing.numberofdofs;i++){ - if(indexing.doftype[i]==approximation_enum){ - outdoflist[count]=indexing.doflist[i]; - count++; - } - } - ISSMASSERT(count); + + if(approximation_enum==NoneApproximationEnum){ + if(setenum==GsetEnum)for(i=0;iindexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i]; + if(setenum==FsetEnum)for(i=0;iindexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i]; + if(setenum==SsetEnum)for(i=0;iindexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i]; } else{ - for(i=0;iindexing.numberofdofs;i++){ - outdoflist[i]=indexing.doflist[i]; - } + + if(setenum==GsetEnum){ + if(indexing.doftype){ + count=0; + for(i=0;iindexing.gsize;i++){ + if(indexing.doftype[i]==approximation_enum){ + outdoflist[count]=indexing.gdoflist[i]; + count++; + } + } + ISSMASSERT(count); + } + else for(i=0;iindexing.gsize;i++) outdoflist[i]=indexing.gdoflist[i]; + } + else if(setenum==FsetEnum){ + if(indexing.doftype){ + count=0; + for(i=0;iindexing.gsize;i++){ + if((indexing.doftype[i]==approximation_enum) && (indexing.f_set[i])){ + outdoflist[count]=indexing.fdoflist[i]; + count++; + } + } + ISSMASSERT(count); + } + else for(i=0;iindexing.fsize;i++) outdoflist[i]=indexing.fdoflist[i]; + } + else if(setenum==SsetEnum){ + if(indexing.doftype){ + count=0; + for(i=0;iindexing.gsize;i++){ + if((indexing.doftype[i]==approximation_enum) && (indexing.s_set[i])){ + outdoflist[count]=indexing.sdoflist[i]; + count++; + } + } + ISSMASSERT(count); + } + else for(i=0;iindexing.ssize;i++) outdoflist[i]=indexing.sdoflist[i]; + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + } +} +/*}}}*/ +/*FUNCTION Node::GetInternalDofList{{{1*/ +void Node::GetInternalDofList(int* outdoflist,int approximation_enum,int setenum){ + int i; + int count=0; + int count2=0; + + if(approximation_enum==NoneApproximationEnum){ + if(setenum==GsetEnum)for(i=0;iindexing.gsize;i++) outdoflist[i]=i; + else if(setenum==FsetEnum){ + count=0; + for(i=0;iindexing.fsize;i++){ + if(indexing.f_set[i]){ + outdoflist[count]=i; + count++; + } + } + } + else if(setenum==SsetEnum){ + count=0; + for(i=0;iindexing.ssize;i++){ + if(indexing.s_set[i]){ + outdoflist[count]=i; + count++; + } + } + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + } + else{ + + if(setenum==GsetEnum){ + if(indexing.doftype){ + count=0; + for(i=0;iindexing.gsize;i++){ + if(indexing.doftype[i]==approximation_enum){ + outdoflist[count]=count; + count++; + } + } + ISSMASSERT(count); + } + else for(i=0;iindexing.gsize;i++) outdoflist[i]=i; + } + else if(setenum==FsetEnum){ + + if(indexing.doftype){ + count=0; + count2=0; + for(i=0;iindexing.gsize;i++){ + if((indexing.doftype[i]==approximation_enum) && (indexing.f_set[i])){ + outdoflist[count]=count2; + count++; + } + if(indexing.doftype[i]==approximation_enum)count2++; + } + ISSMASSERT(count); + } + else{ + + count=0; + for(i=0;iindexing.gsize;i++){ + if(indexing.f_set[i]){ + outdoflist[count]=i; + count++; + } + } + } + } + else if(setenum==SsetEnum){ + if(indexing.doftype){ + count=0; + count2=0; + for(i=0;iindexing.gsize;i++){ + if((indexing.doftype[i]==approximation_enum) && (indexing.s_set[i])){ + outdoflist[count]=count2; + count++; + } + if(indexing.doftype[i]==approximation_enum)count2++; + } + ISSMASSERT(count); + } + else{ + count=0; + for(i=0;iindexing.gsize;i++){ + if(indexing.s_set[i]){ + outdoflist[count]=i; + count++; + } + } + } + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); } } @@ -380,24 +518,12 @@ /*Node numerics:*/ /*FUNCTION Node::ApplyConstraints{{{1*/ -void Node::ApplyConstraint(Vec yg,int dof,double value){ +void Node::ApplyConstraint(int dof,double value){ int index; - /*First, dof should be added in the s set, describing which + /*Dof should be added in the s set, describing which * dofs are constrained to a certain value (dirichlet boundary condition*/ - DofInSSet(dof-1); - - /*Second, we should add value into yg, at dof corresponding to doflist[dof], unless - * we are a clone!*/ - - if(!indexing.clone){ - - index=indexing.doflist[dof-1]; //matlab indexing - - VecSetValues(yg,1,&index,&value,INSERT_VALUES); - - } - + this->indexing.svalues[dof-1]=value; } /*}}}*/ @@ -410,18 +536,44 @@ int i; - for(i=0;iindexing.numberofdofs;i++){ + for(i=0;iindexing.gsize;i++){ /*g set: */ - VecSetValues(pv_g,1,&indexing.doflist[i],&gvalue,INSERT_VALUES); + VecSetValues(pv_g,1,&indexing.gdoflist[i],&gvalue,INSERT_VALUES); /*f set: */ value=(double)this->indexing.f_set[i]; - VecSetValues(pv_f,1,&indexing.doflist[i],&value,INSERT_VALUES); + VecSetValues(pv_f,1,&indexing.gdoflist[i],&value,INSERT_VALUES); /*s set: */ value=(double)this->indexing.s_set[i]; - VecSetValues(pv_s,1,&indexing.doflist[i],&value,INSERT_VALUES); - - } + VecSetValues(pv_s,1,&indexing.gdoflist[i],&value,INSERT_VALUES); + + } + + +} +/*}}}*/ +/*FUNCTION Node::CreateNodalConstraints{{{1*/ +void Node::CreateNodalConstraints(Vec ys){ + + int i; + double* values=NULL; + int count; + + /*Recover values for s set and plug them in constraints vector: */ + if(this->indexing.ssize){ + values=(double*)xmalloc(this->indexing.ssize*sizeof(double)); + count=0; + for(i=0;iindexing.gsize;i++){ + if(this->indexing.s_set[i])values[count]=this->indexing.svalues[i]; + count++; + } + + /*Add values into constraint vector: */ + VecSetValues(ys,this->indexing.ssize,this->indexing.sdoflist,values,INSERT_VALUES); + } + + /*Free ressources:*/ + xfree((void**)&values); @@ -467,20 +619,49 @@ /*}}}*/ /*FUNCTION Node::GetNumberOfDofs{{{1*/ -int Node::GetNumberOfDofs(int approximation_enum){ +int Node::GetNumberOfDofs(int approximation_enum,int setenum){ + + /*Get number of degrees of freedom in a node, for a certain set (g,f or s-set) + *and for a certain approximation type: */ int i; int numdofs=0; - /*Count the dofs if an approximation is specified and the node contains several dofs type*/ - if(approximation_enum && this->indexing.doftype){ - for(i=0;iindexing.numberofdofs;i++){ - if(this->indexing.doftype[i]==approximation_enum) numdofs++; - } - } - else numdofs=this->indexing.numberofdofs; - - ISSMASSERT(numdofs); + if(approximation_enum==NoneApproximationEnum){ + if (setenum==GsetEnum) numdofs=this->indexing.gsize; + else if (setenum==FsetEnum) numdofs=this->indexing.fsize; + else if (setenum==SsetEnum) numdofs=this->indexing.ssize; + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + } + else{ + if(setenum==GsetEnum){ + if(this->indexing.doftype){ + numdofs=0; + for(i=0;iindexing.gsize;i++){ + if(this->indexing.doftype[i]==approximation_enum) numdofs++; + } + } + else numdofs=this->indexing.gsize; + } + else if (setenum==FsetEnum){ + if(this->indexing.doftype){ + numdofs=0; + for(i=0;iindexing.gsize;i++){ + if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.f_set[i])) numdofs++; + } + } + else numdofs=this->indexing.fsize; + } + else if (setenum==SsetEnum){ + if(this->indexing.doftype){ + numdofs=0; + for(i=0;iindexing.gsize;i++){ + if((this->indexing.doftype[i]==approximation_enum) && (this->indexing.s_set[i])) numdofs++; + } + } + else numdofs=this->indexing.ssize; + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + } return numdofs; - } /*}}}*/ @@ -625,5 +806,5 @@ /* DofObject routines:*/ /*FUNCTION Node::DistributeDofs{{{1*/ -void Node::DistributeDofs(int* pdofcount){ +void Node::DistributeDofs(int* pdofcount,int setenum){ int i; @@ -638,9 +819,28 @@ } - /*This node should distribute dofs, go ahead: */ - for(i=0;iindexing.numberofdofs;i++){ - indexing.doflist[i]=dofcount+i; - } - dofcount+=this->indexing.numberofdofs; + /*This node should distribute dofs for setenum set (eg, f_set or s_set), go ahead: */ + + if(setenum==GsetEnum){ + for(i=0;iindexing.gsize;i++){ + indexing.gdoflist[i]=dofcount+i; + } + dofcount+=this->indexing.gsize; + } + else if(setenum==FsetEnum){ + this->indexing.InitSet(setenum); + for(i=0;iindexing.fsize;i++){ + indexing.fdoflist[i]=dofcount+i; + } + dofcount+=this->indexing.fsize; + } + else if(setenum==SsetEnum){ + this->indexing.InitSet(setenum); + for(i=0;iindexing.ssize;i++){ + indexing.sdoflist[i]=dofcount+i; + } + dofcount+=this->indexing.ssize; + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + /*Assign output pointers: */ @@ -649,6 +849,6 @@ } /*}}}*/ -/*FUNCTION Node::OffsetDofs{{{1*/ -void Node::OffsetDofs(int dofcount){ +/*FUNCTION Node::Off_setDofs{{{1*/ +void Node::OffsetDofs(int dofcount,int setenum){ int i; @@ -656,16 +856,23 @@ if(indexing.clone){ - /*This node is a clone, don't offset the dofs!: */ + /*This node is a clone, don't off_set the dofs!: */ return; } - /*This node should offset the dofs, go ahead: */ - for(i=0;iindexing.numberofdofs;i++){ - indexing.doflist[i]+=dofcount; - } + /*This node should off_set the dofs, go ahead: */ + if(setenum==GsetEnum){ + for(i=0;iindexing.gsize;i++) indexing.gdoflist[i]+=dofcount; + } + else if(setenum==FsetEnum){ + for(i=0;iindexing.fsize;i++) indexing.fdoflist[i]+=dofcount; + } + else if(setenum==SsetEnum){ + for(i=0;iindexing.ssize;i++) indexing.sdoflist[i]+=dofcount; + } + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); } /*}}}*/ /*FUNCTION Node::ShowTrueDofs{{{1*/ -void Node::ShowTrueDofs(int* truedofs, int ncols){ +void Node::ShowTrueDofs(int* truedofs, int ncols,int setenum){ int j; @@ -676,11 +883,13 @@ /*Ok, we are not a clone, just plug our dofs into truedofs: */ - for(j=0;jindexing.numberofdofs;j++){ - *(truedofs+ncols*sid+j)=indexing.doflist[j]; - } + if(setenum==GsetEnum)for(j=0;jindexing.gsize;j++) *(truedofs+ncols*sid+j)=indexing.gdoflist[j]; + else if(setenum==FsetEnum)for(j=0;jindexing.fsize;j++) *(truedofs+ncols*sid+j)=indexing.fdoflist[j]; + else if(setenum==SsetEnum)for(j=0;jindexing.ssize;j++) *(truedofs+ncols*sid+j)=indexing.sdoflist[j]; + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); + } /*}}}*/ /*FUNCTION Node::UpdateCloneDofs{{{1*/ -void Node::UpdateCloneDofs(int* alltruedofs,int ncols){ +void Node::UpdateCloneDofs(int* alltruedofs,int ncols,int setenum){ int j; @@ -693,7 +902,8 @@ /*Ok, we are a clone node, but we did not create the dofs for this node. * * Therefore, our doflist is garbage right now. Go pick it up in the alltruedofs: */ - for(j=0;jindexing.numberofdofs;j++){ - indexing.doflist[j]=*(alltruedofs+ncols*sid+j); - } + if(setenum==GsetEnum)for(j=0;jindexing.gsize;j++) indexing.gdoflist[j]=*(alltruedofs+ncols*sid+j); + else if(setenum==FsetEnum)for(j=0;jindexing.fsize;j++) indexing.fdoflist[j]=*(alltruedofs+ncols*sid+j); + else if(setenum==SsetEnum)for(j=0;jindexing.ssize;j++) indexing.sdoflist[j]=*(alltruedofs+ncols*sid+j); + else ISSMERROR("%s%s%s"," set of enum type ",EnumToString(setenum)," not supported yet!"); } Index: /issm/trunk/src/c/objects/Node.h =================================================================== --- /issm/trunk/src/c/objects/Node.h (revision 5771) +++ /issm/trunk/src/c/objects/Node.h (revision 5772) @@ -65,4 +65,5 @@ /*Node numerical routines {{{1*/ void Configure(DataSet* nodes,Vertices* vertices); + void CreateNodalConstraints(Vec ys); void SetCurrentConfiguration(DataSet* nodes,Vertices* vertices); int Sid(void); @@ -72,12 +73,13 @@ bool InAnalysis(int analysis_type); int GetApproximation(); - int GetNumberOfDofs(int approximation_enum=0); + int GetNumberOfDofs(int approximation_enum,int setenum); int IsClone(); - void ApplyConstraint(Vec yg,int dof,double value); + void ApplyConstraint(int dof,double value); void DofInSSet(int dof); - int GetDof(int dofindex); + int GetDof(int dofindex,int setenum); void CreateVecSets(Vec pv_g,Vec pv_f,Vec pv_s); int GetConnectivity(); - void GetDofList(int* poutdoflist,int approximation_enum=0); + void GetDofList(int* poutdoflist,int approximation_enum,int setenum); + void GetInternalDofList(int* poutdoflist,int approximation_enum,int setenum); int GetDofList1(void); double GetX(); @@ -92,8 +94,8 @@ /*}}}*/ /*Dof Object routines {{{1*/ - void DistributeDofs(int* pdofcount); - void OffsetDofs(int dofcount); - void ShowTrueDofs(int* truerows,int ncols); - void UpdateCloneDofs(int* alltruerows,int ncols); + void DistributeDofs(int* pdofcount,int setenum); + void OffsetDofs(int dofcount,int setenum); + void ShowTrueDofs(int* truerows,int ncols,int setenum); + void UpdateCloneDofs(int* alltruerows,int ncols,int setenum); void SetClone(int* minranks); void CreatePartition(Vec partition); Index: /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5772) +++ /issm/trunk/src/c/objects/Numerics/ElementMatrix.cpp (revision 5772) @@ -0,0 +1,263 @@ +/*!\file ElementMatrix.cpp + * \brief: implementation of the ElementMatrix object, used to plug values from element into global stiffness matrix + */ + +/*Headers:*/ +/*{{{1*/ +#ifdef HAVE_CONFIG_H + #include "config.h" +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "stdio.h" +#include +#include "../objects.h" +#include "../../shared/shared.h" +#include "../../Container/Container.h" +#include "../../include/include.h" +/*}}}*/ + +/*ElementMatrix constructors and destructor*/ +/*FUNCTION ElementMatrix::ElementMatrix(){{{1*/ +ElementMatrix::ElementMatrix(){ + + this->nrows=0; + this->ncols=0; + this->symmetric=false; + this->row_fsize=0; + this->values=NULL; + this->row_finternaldoflist=NULL; + this->row_fexternaldoflist=NULL; + this->row_ssize=0; + this->row_sinternaldoflist=NULL; + this->row_sexternaldoflist=NULL; + this->col_fsize=0; + this->col_finternaldoflist=NULL; + this->col_fexternaldoflist=NULL; + this->col_ssize=0; + this->col_sinternaldoflist=NULL; + this->col_sexternaldoflist=NULL; + +} +/*}}}*/ +/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* in_gexternaldoflist){{{1*/ +ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* in_gexternaldoflist){ + + if(symmetric=false)ISSMERROR(" calling symmetric constructor with false symmetric flag!"); + + this->nrows=gsize; + this->ncols=gsize; + this->symmetric=true; + this->kff=false; + + /*fill values with 0: */ + this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double)); + + if(this->nrows){ + this->gexternaldoflist=(int*)xmalloc(this->nrows*sizeof(int)); + memcpy(this->gexternaldoflist,in_gexternaldoflist,this->nrows*sizeof(int)); + } + else this->gexternaldoflist=NULL; + + /*don't do rows and cols, because kff is false:*/ + this->row_fsize=0; + this->row_finternaldoflist=NULL; + this->row_fexternaldoflist=NULL; + this->row_ssize=0; + this->row_sinternaldoflist=NULL; + this->row_sexternaldoflist=NULL; + + this->col_fsize=0; + this->col_finternaldoflist=NULL; + this->col_fexternaldoflist=NULL; + this->col_ssize=0; + this->col_sinternaldoflist=NULL; + this->col_sexternaldoflist=NULL; + +} +/*}}}*/ +/*FUNCTION ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){{{1*/ +ElementMatrix::ElementMatrix(int gsize,bool symmetric,int* finternaldoflist,int* fexternaldoflist,int fsize,int* sinternaldoflist,int* sexternaldoflist,int ssize){ + + if(symmetric=false)ISSMERROR(" calling symmetric constructor with false symmetric flag!"); + + this->nrows=gsize; + this->ncols=gsize; + this->symmetric=true; + this->kff=false; + + /*fill values with 0: */ + this->values=(double*)xcalloc(this->nrows*this->ncols,sizeof(double)); + + /*row dofs: */ + this->row_fsize=fsize; + if(fsize){ + this->row_finternaldoflist=(int*)xmalloc(fsize*sizeof(int)); + this->row_fexternaldoflist=(int*)xmalloc(fsize*sizeof(int)); + memcpy(this->row_finternaldoflist,finternaldoflist,fsize*sizeof(int)); + memcpy(this->row_fexternaldoflist,fexternaldoflist,fsize*sizeof(int)); + } + else{ + this->row_finternaldoflist=NULL; + this->row_fexternaldoflist=NULL; + } + + this->row_ssize=ssize; + if(ssize){ + this->row_sinternaldoflist=(int*)xmalloc(ssize*sizeof(int)); + this->row_sexternaldoflist=(int*)xmalloc(ssize*sizeof(int)); + memcpy(this->row_sinternaldoflist,sinternaldoflist,ssize*sizeof(int)); + memcpy(this->row_sexternaldoflist,sexternaldoflist,ssize*sizeof(int)); + } + else{ + this->row_sinternaldoflist=NULL; + this->row_sexternaldoflist=NULL; + } + + /*don't do cols, we can't pick them up from the rows: */ + this->col_fsize=0; + this->col_finternaldoflist=NULL; + this->col_fexternaldoflist=NULL; + this->col_ssize=0; + this->col_sinternaldoflist=NULL; + this->col_sexternaldoflist=NULL; + + /*don't do g-set: */ + this->gexternaldoflist=NULL; + +} +/*}}}*/ +/*FUNCTION ElementMatrix::~ElementMatrix(){{{1*/ +ElementMatrix::~ElementMatrix(){ + + xfree((void**)&this->values); + xfree((void**)&this->gexternaldoflist); + xfree((void**)&this->row_finternaldoflist); + xfree((void**)&this->row_fexternaldoflist); + xfree((void**)&this->row_sinternaldoflist); + xfree((void**)&this->row_sexternaldoflist); + xfree((void**)&this->col_finternaldoflist); + xfree((void**)&this->col_fexternaldoflist); + xfree((void**)&this->col_sinternaldoflist); + xfree((void**)&this->col_sexternaldoflist); +} +/*}}}*/ + +/*ElementMatrix specific routines: */ +/*FUNCTION ElementMatrix::AddValues(double* Ke_gg){{{1*/ +void ElementMatrix::AddValues(double* Ke_gg){ + + int i,j; + + if(Ke_gg){ + for (i=0;inrows;i++){ + for(j=0;jncols;j++){ + *(this->values+this->ncols*i+j)+=*(Ke_gg+this->ncols*i+j); + } + } + } +} +/*}}}*/ +/*FUNCTION ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){{{1*/ +void ElementMatrix::AddToGlobal(Mat Kgg, Mat Kff, Mat Kfs){ + + int i,j; + double* internalvalues=NULL; + + if(this->symmetric){ + /*only use row dofs to add values into global matrices: */ + + if(!this->kff){ + /*add internal values into global matrix, using the fexternaldoflist: */ + MatSetValues(Kgg,this->nrows,this->gexternaldoflist,this->nrows,this->gexternaldoflist,(const double*)values,ADD_VALUES); + } + else{ + if(this->row_fsize){ + /*first, retrieve values that are in the f-set from the g-set values matrix: */ + internalvalues=(double*)xmalloc(this->row_fsize*this->row_fsize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_fsize;j++){ + *(internalvalues+this->row_fsize*i+j)=*(this->values+this->ncols*this->row_finternaldoflist[i]+this->row_finternaldoflist[j]); + } + } + /*add internal values into global matrix, using the fexternaldoflist: */ + MatSetValues(Kff,this->row_fsize,this->row_fexternaldoflist,this->row_fsize,this->row_fexternaldoflist,(const double*)internalvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&internalvalues); + } + + + if(this->row_ssize){ + /*first, retrieve values that are in the f and s-set from the g-set values matrix: */ + internalvalues=(double*)xmalloc(this->row_fsize*this->row_ssize*sizeof(double)); + for(i=0;irow_fsize;i++){ + for(j=0;jrow_ssize;j++){ + *(internalvalues+this->row_ssize*i+j)=*(this->values+this->ncols*this->row_finternaldoflist[i]+this->row_sinternaldoflist[j]); + } + } + /*add internal values into global matrix, using the fexternaldoflist: */ + MatSetValues(Kfs,this->row_fsize,this->row_fexternaldoflist,this->row_ssize,this->row_sexternaldoflist,(const double*)internalvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&internalvalues); + } + } + } + else{ + ISSMERROR(" unsymmetric AddToGlobal routine not support yet!"); + } + +} +/*}}}*/ +/*FUNCTION ElementMatrix::Echo(void){{{1*/ +void ElementMatrix::Echo(void){ + + int i,j; + printf("Element Matrix echo: \n"); + printf(" nrows: %i\n",nrows); + printf(" ncols: %i\n",ncols); + printf(" symmetric: %s\n",symmetric?"true":"false"); + printf(" kff: %s\n",kff?"true":"false"); + + printf(" values: \n"); + for(i=0;i +#include "../objects.h" +#include "../../shared/shared.h" +#include "../../Container/Container.h" +#include "../../include/include.h" +/*}}}*/ + +/*ElementVector constructors and destructor*/ +/*FUNCTION ElementVector::ElementVector(){{{1*/ +ElementVector::ElementVector(){ + + this->nrows=0; + this->values=NULL; + this->fsize=0; + this->finternaldoflist=NULL; + this->fexternaldoflist=NULL; + +} +/*}}}*/ +/*FUNCTION ElementVector::ElementVector(int gsize,int* gexternaldoflist){{{1*/ +ElementVector::ElementVector(int gsize,int* in_gexternaldoflist){ + + this->nrows=gsize; + this->pf=false; + + /*fill values with 0: */ + this->values=(double*)xcalloc(this->nrows,sizeof(double)); + + /*dofs: */ + if(this->nrows){ + this->gexternaldoflist=(int*)xmalloc(nrows*sizeof(int)); + memcpy(this->gexternaldoflist,in_gexternaldoflist,nrows*sizeof(int)); + } + else{ + this->gexternaldoflist=NULL; + } + /*not needed: */ + this->finternaldoflist=NULL; + this->fexternaldoflist=NULL; + +} +/*}}}*/ +/*FUNCTION ElementVector::ElementVector(int gsize,int* finternaldoflist,int* fexternaldoflist,int fsize){{{1*/ +ElementVector::ElementVector(int gsize,int* in_finternaldoflist,int* in_fexternaldoflist,int in_fsize){ + + this->nrows=gsize; + this->pf=true; + + /*fill values with 0: */ + this->values=(double*)xcalloc(this->nrows,sizeof(double)); + + /*dofs: */ + this->fsize=in_fsize; + if(this->fsize){ + this->finternaldoflist=(int*)xmalloc(fsize*sizeof(int)); + this->fexternaldoflist=(int*)xmalloc(fsize*sizeof(int)); + memcpy(this->finternaldoflist,in_finternaldoflist,fsize*sizeof(int)); + memcpy(this->fexternaldoflist,in_fexternaldoflist,fsize*sizeof(int)); + } + else{ + this->finternaldoflist=NULL; + this->fexternaldoflist=NULL; + } + + /*not needed: */ + this->gexternaldoflist=NULL; + +} +/*}}}*/ +/*FUNCTION ElementVector::~ElementVector(){{{1*/ +ElementVector::~ElementVector(){ + + xfree((void**)&this->values); + xfree((void**)&this->gexternaldoflist); + xfree((void**)&this->finternaldoflist); + xfree((void**)&this->fexternaldoflist); +} +/*}}}*/ + +/*ElementVector specific routines: */ +/*FUNCTION ElementVector::AddValues(double* Ke_gg){{{1*/ +void ElementVector::AddValues(double* pe_gg){ + + int i,j; + + if(pe_gg){ + for (i=0;inrows;i++){ + this->values[i]+=pe_gg[i]; + } + } +} +/*}}}*/ +/*FUNCTION ElementVector::AddToGlobal(Vec pg, Vec pf){{{1*/ +void ElementVector::AddToGlobal(Vec pg, Vec pf){ + + int i; + double* internalvalues=NULL; + + if(!pf){ + /*add internal values into global vector, using the fexternaldoflist: */ + VecSetValues(pg,this->nrows,this->gexternaldoflist,(const double*)values,ADD_VALUES); + } + else{ + if(this->fsize){ + /*first, retrieve values that are in the f-set from the g-set values vector: */ + internalvalues=(double*)xmalloc(this->fsize*sizeof(double)); + for(i=0;ifsize;i++){ + internalvalues[i]=this->values[this->finternaldoflist[i]]; + } + /*add internal values into global vector, using the fexternaldoflist: */ + VecSetValues(pf,this->fsize,this->fexternaldoflist,(const double*)internalvalues,ADD_VALUES); + + /*Free ressources:*/ + xfree((void**)&internalvalues); + } + } +} +/*}}}*/ Index: /issm/trunk/src/c/objects/Numerics/ElementVector.h =================================================================== --- /issm/trunk/src/c/objects/Numerics/ElementVector.h (revision 5772) +++ /issm/trunk/src/c/objects/Numerics/ElementVector.h (revision 5772) @@ -0,0 +1,45 @@ +/*!\file: ElementVector.h + * \brief container for information needed to plug element vector generated by elements + * into the pf global load vector. + * This object will hold the element vector on the g-set, the internal as well as external + * dof lists in the f set + */ + +#ifndef _ELEMENT_VECTOR_H_ +#define _ELEMENT_VECTOR_H_ + +/*Headers:*/ +/*{{{1*/ +#include "../Object.h" +#include "../../toolkits/toolkits.h" +/*}}}*/ + +class ElementVector{ + + public: + + int nrows; + double* values; + bool pf; + + //gset + int* gexternaldoflist; + + //fset + int fsize; + int* finternaldoflist; + int* fexternaldoflist; + + /*ElementVector constructors, destructors {{{1*/ + ElementVector(); + ElementVector(int gsize,int* gexternaldoflist); + ElementVector(int gsize,int* finternaldoflist,int* fexternaldoflist,int fsize); + ~ElementVector(); + /*}}}*/ + /*ElementVector specific routines {{{1*/ + void AddValues(double* pe_gg); + void AddToGlobal(Vec pg, Vec pf); + /*}}}*/ +}; +#endif //#ifndef _ELEMENT_VECTOR_H_ + Index: /issm/trunk/src/c/objects/objects.h =================================================================== --- /issm/trunk/src/c/objects/objects.h (revision 5771) +++ /issm/trunk/src/c/objects/objects.h (revision 5772) @@ -73,4 +73,8 @@ #include "./Materials/Matpar.h" +/*Numerics:*/ +#include "./Numerics/ElementMatrix.h" +#include "./Numerics/ElementVector.h" + /*Params: */ Index: /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 5771) +++ /issm/trunk/src/c/solutions/ResetBoundaryConditions.cpp (revision 5772) @@ -10,5 +10,6 @@ int verbose=0; - Vec ug=NULL; + Vec yg=NULL; + Vec ys=NULL; int analysis_counter; @@ -19,22 +20,17 @@ femmodel->SetCurrentConfiguration(analysis_type); - GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + GetSolutionFromInputsx( &yg, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); - /*For this analysis_type, free existing boundary condition vectors: */ + /*For this analysis_type, free existing boundary condition vector: */ analysis_counter=femmodel->analysis_counter; - - //global dof set - VecFree(&femmodel->m_yg[analysis_counter]); - //in the s-set VecFree(&femmodel->m_ys[analysis_counter]); - //Now, duplicate ug (the solution vector) into the boundary conditions vector on the g-set - VecDuplicatePatch(&femmodel->m_yg[analysis_counter],ug); - //Reduce from g to s set - Reducevectorgtosx(&femmodel->m_ys[analysis_counter],femmodel->m_yg[analysis_counter],femmodel->m_nodesets[analysis_counter]); + Reducevectorgtosx(&ys,yg,femmodel->m_nodesets[analysis_counter]); + + /*Plug into femmodel->m_ys: */ + femmodel->m_ys[analysis_counter]=ys; /*Free ressources:*/ - VecFree(&ug); - + VecFree(&yg); } Index: /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp (revision 5771) +++ /issm/trunk/src/c/solvers/solver_adjoint_linear.cpp (revision 5772) @@ -17,5 +17,5 @@ Vec pg = NULL, pf=NULL; - SystemMatricesx(&Kgg, &pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 5771) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp (revision 5772) @@ -12,90 +12,12 @@ void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads){ - /*intermediary: */ - Mat Kgg = NULL, Kff = NULL, Kfs = NULL; - Vec ug = NULL, uf = NULL, old_ug= NULL, old_uf = NULL; - Vec pg = NULL, pf = NULL; - Loads* loads=NULL; - int converged; - int constraints_converged; - int num_unstable_constraints; - int count; - int min_mechanical_constraints; - int max_nonlinear_iterations; - /*parameters:*/ - int verbose=0; + bool kff=false; /*Recover parameters: */ - femmodel->parameters->FindParam(&verbose,VerboseEnum); - femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); - femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); - - /*Were loads requested as output? : */ - if(conserve_loads){ - loads=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution - } - else{ - loads=(Loads*)femmodel->loads; //modify loads in this solution - } + femmodel->parameters->FindParam(&kff,KffEnum); - count=1; - converged=0; + if(!kff) solver_diagnostic_nonlinear_kgg(femmodel,conserve_loads); + else solver_diagnostic_nonlinear_kff(femmodel,conserve_loads); - /*Start non-linear iteration using input velocity: */ - GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); - Reducevectorgtofx(&uf, ug, femmodel->nodesets); - - //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) - InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - - for(;;){ - - //save pointer to old velocity - VecFree(&old_ug);old_ug=ug; - VecFree(&old_uf);old_uf=uf; - - SystemMatricesx(&Kgg, &pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); - - Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); - - Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg); MatFree(&Kfs); - - Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); - - Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters); - - InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); - - PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); - if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); - - convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); - - InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); - - //rift convergence - if (!constraints_converged) { - if (converged){ - if (num_unstable_constraints <= min_mechanical_constraints) converged=1; - else converged=0; - } - } - - /*Increase count: */ - count++; - if(converged==1)break; - if(count>=max_nonlinear_iterations){ - _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); - break; - } - } - - /*clean-up*/ - if(conserve_loads) delete loads; - VecFree(&uf); - VecFree(&ug); - VecFree(&old_uf); - VecFree(&old_ug); - } Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp (revision 5772) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kff.cpp (revision 5772) @@ -0,0 +1,99 @@ +/*!\file: solver_diagnostic_nonlinear_kff.cpp + * \brief: core of the diagnostic solution for non linear materials, using f and s-set + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" +#include "../solutions/solutions.h" +#include "./solvers.h" + +void solver_diagnostic_nonlinear_kff(FemModel* femmodel,bool conserve_loads){ + + /*intermediary: */ + Mat Kff = NULL, Kfs = NULL; + Vec ug = NULL, uf = NULL, old_ug= NULL, old_uf = NULL; + Vec pf = NULL; + Loads* loads=NULL; + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int min_mechanical_constraints; + int max_nonlinear_iterations; + + /*parameters:*/ + int verbose=0; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + + /*Were loads requested as output? : */ + if(conserve_loads){ + loads=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution + } + else{ + loads=(Loads*)femmodel->loads; //modify loads in this solution + } + + count=1; + converged=0; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf, ug, femmodel->nodesets); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + for(;;){ + + //save pointer to old velocity + VecFree(&old_ug);old_ug=ug; + VecFree(&old_uf);old_uf=uf; + + SystemMatricesx(NULL,&Kff, &Kfs, NULL,&pf, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + + Reduceloadx(pf, Kfs, femmodel->ys, femmodel->parameters); MatFree(&Kfs); + + Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); + + Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters); + + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); + + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); + + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + + //rift convergence + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=1; + else converged=0; + } + } + + /*Increase count: */ + count++; + if(converged==1)break; + if(count>=max_nonlinear_iterations){ + _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + break; + } + } + + /*clean-up*/ + if(conserve_loads) delete loads; + VecFree(&uf); + VecFree(&ug); + VecFree(&old_uf); + VecFree(&old_ug); + +} Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp (revision 5772) +++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear_kgg.cpp (revision 5772) @@ -0,0 +1,101 @@ +/*!\file: solver_diagnostic_nonlinear.cpp + * \brief: core of the diagnostic solution for non linear materials + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "../modules/modules.h" +#include "../solutions/solutions.h" +#include "./solvers.h" + +void solver_diagnostic_nonlinear_kgg(FemModel* femmodel,bool conserve_loads){ + + /*intermediary: */ + Mat Kgg = NULL, Kff = NULL, Kfs = NULL; + Vec ug = NULL, uf = NULL, old_ug= NULL, old_uf = NULL; + Vec pg = NULL, pf = NULL; + Loads* loads=NULL; + int converged; + int constraints_converged; + int num_unstable_constraints; + int count; + int min_mechanical_constraints; + int max_nonlinear_iterations; + + /*parameters:*/ + int verbose=0; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&verbose,VerboseEnum); + femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum); + femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum); + + /*Were loads requested as output? : */ + if(conserve_loads){ + loads=(Loads*)femmodel->loads->Copy(); //protect loads from being modified by the solution + } + else{ + loads=(Loads*)femmodel->loads; //modify loads in this solution + } + + count=1; + converged=0; + + /*Start non-linear iteration using input velocity: */ + GetSolutionFromInputsx(&ug, femmodel->elements, femmodel->nodes, femmodel->vertices, loads, femmodel->materials, femmodel->parameters); + Reducevectorgtofx(&uf, ug, femmodel->nodesets); + + //Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromSolutionx(femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + for(;;){ + + //save pointer to old velocity + VecFree(&old_ug);old_ug=ug; + VecFree(&old_uf);old_uf=uf; + + SystemMatricesx(&Kgg, NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + + Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); + + Reduceloadfromgtofx(&pf, pg, Kfs, femmodel->ys, femmodel->nodesets,femmodel->parameters); VecFree(&pg); MatFree(&Kfs); + + Solverx(&uf, Kff, pf, old_uf, femmodel->parameters); + + Mergesolutionfromftogx(&ug, uf,femmodel->ys,femmodel->nodesets,femmodel->parameters); + + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug); + + PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, femmodel->elements,femmodel->nodes,femmodel->vertices,loads,femmodel->materials,femmodel->parameters); + if(verbose)_printf_(" number of unstable constraints: %i\n",num_unstable_constraints); + + convergence(&converged,Kff,pf,uf,old_uf,femmodel->parameters); MatFree(&Kff);VecFree(&pf); + + InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum); + + //rift convergence + if (!constraints_converged) { + if (converged){ + if (num_unstable_constraints <= min_mechanical_constraints) converged=1; + else converged=0; + } + } + + /*Increase count: */ + count++; + if(converged==1)break; + if(count>=max_nonlinear_iterations){ + _printf_(" maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations); + break; + } + } + + /*clean-up*/ + if(conserve_loads) delete loads; + VecFree(&uf); + VecFree(&ug); + VecFree(&old_uf); + VecFree(&old_ug); + +} Index: /issm/trunk/src/c/solvers/solver_linear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_linear.cpp (revision 5771) +++ /issm/trunk/src/c/solvers/solver_linear.cpp (revision 5772) @@ -15,5 +15,5 @@ Vec pg = NULL, pf = NULL; - SystemMatricesx(&Kgg,&pg,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); Index: /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp =================================================================== --- /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 5771) +++ /issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp (revision 5772) @@ -51,5 +51,5 @@ - SystemMatricesx(&Kgg,&pg,&melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + SystemMatricesx(&Kgg,NULL, NULL, &pg,NULL, &melting_offset,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); Reducematrixfromgtofx(&Kff,&Kfs,Kgg,femmodel->nodesets,femmodel->parameters); MatFree(&Kgg); Index: /issm/trunk/src/c/solvers/solvers.h =================================================================== --- /issm/trunk/src/c/solvers/solvers.h (revision 5771) +++ /issm/trunk/src/c/solvers/solvers.h (revision 5772) @@ -14,4 +14,6 @@ void solver_thermal_nonlinear(FemModel* femmodel); void solver_diagnostic_nonlinear(FemModel* femmodel,bool conserve_loads); +void solver_diagnostic_nonlinear_kgg(FemModel* femmodel,bool conserve_loads); +void solver_diagnostic_nonlinear_kff(FemModel* femmodel,bool conserve_loads); void solver_linear(FemModel* femmodel); void solver_adjoint_linear(FemModel* femmodel); Index: /issm/trunk/src/m/classes/@model/model.m =================================================================== --- /issm/trunk/src/m/classes/@model/model.m (revision 5771) +++ /issm/trunk/src/m/classes/@model/model.m (revision 5772) @@ -309,4 +309,7 @@ md.qmu_relax=0; + %solution parameters + md.kff=0; + %partitioner: md.adjacency=[]; Index: /issm/trunk/src/m/classes/@model/setdefaultparameters.m =================================================================== --- /issm/trunk/src/m/classes/@model/setdefaultparameters.m (revision 5771) +++ /issm/trunk/src/m/classes/@model/setdefaultparameters.m (revision 5772) @@ -256,2 +256,5 @@ %Ice solver: 'general' for Matlab's default solver (or 'lu' or 'sholesky') md.solver_type='general'; + +%solution speed-up +md.kff=0; Index: /issm/trunk/src/m/classes/public/marshall.m =================================================================== --- /issm/trunk/src/m/classes/public/marshall.m (revision 5771) +++ /issm/trunk/src/m/classes/public/marshall.m (revision 5772) @@ -138,4 +138,5 @@ WriteData(fid,md.stokesreconditioning,'Scalar','stokesreconditioning'); WriteData(fid,md.waitonlock,'Scalar','waitonlock'); +WriteData(fid,md.kff,'Integer','kff'); %Thermal parameters