Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 15688) @@ -3268,5 +3268,5 @@ penta_node_ids[14]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[9*index+8]+1; break; - case P1P1Enum: case P1P1GLSEnum: case MINIcondensedEnum: + case P1P1Enum: case P1P1GLSEnum: numnodes = 12; penta_node_ids = xNew(numnodes); @@ -3285,5 +3285,5 @@ penta_node_ids[11]=iomodel->nodecounter+iomodel->numberofvertices+iomodel->elements[6*index+5]; break; - case MINIEnum: + case MINIEnum: case MINIcondensedEnum: numnodes = 13; penta_node_ids = xNew(numnodes); @@ -5792,5 +5792,5 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocityFinal(); + int vnumnodes = this->NumberofNodesVelocity(); int pnumnodes = this->NumberofNodesPressure(); int vnumdof = vnumnodes*NDOF3; @@ -7418,4 +7418,10 @@ Ke =new ElementMatrix(Ke1,Ke2); + /*Condense if requested*/ + if(this->element_type==MINIcondensedEnum){ + int indices[3]={18,19,20}; + Ke->StaticCondensation(3,&indices[0]); + } + /*clean-up and return*/ delete Ke1; @@ -8441,4 +8447,16 @@ pe =new ElementVector(pe1,pe2,pe3); + /*Condense if requested*/ + if(this->element_type==MINIcondensedEnum){ + int indices[3]={18,19,20}; + + this->element_type=MINIEnum; + ElementMatrix* Ke = CreateKMatrixDiagnosticFS(); + this->element_type=MINIcondensedEnum; + + pe->StaticCondensation(Ke,3,&indices[0]); + delete Ke; + } + /*clean-up and return*/ delete pe1; @@ -9250,5 +9268,5 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocityFinal(); + int vnumnodes = this->NumberofNodesVelocity(); int pnumnodes = this->NumberofNodesPressure(); int vnumdof = vnumnodes*NDOF3; @@ -10099,5 +10117,5 @@ /*Fetch number of nodes and dof for this finite element*/ - int vnumnodes = this->NumberofNodesVelocityFinal(); + int vnumnodes = this->NumberofNodesVelocity(); int pnumnodes = this->NumberofNodesPressure(); int vnumdof = vnumnodes*NDOF3; Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp (revision 15688) @@ -1929,5 +1929,5 @@ case P1P1Enum: return NUMNODESP1*2; case P1P1GLSEnum: return NUMNODESP1*2; - case MINIcondensedEnum: return NUMNODESP1*2; + case MINIcondensedEnum: return NUMNODESP1b+NUMNODESP1; case MINIEnum: return NUMNODESP1b+NUMNODESP1; case TaylorHoodEnum: return NUMNODESP2+NUMNODESP1; @@ -1968,19 +1968,2 @@ } /*}}}*/ -/*FUNCTION PentaRef::NumberofNodesVelocityFinal{{{*/ -int PentaRef::NumberofNodesVelocityFinal(void){ - - /*When static condensation is applied, the final number of nodes might be different*/ - - switch(this->element_type){ - case P1P1Enum: return NUMNODESP1; - case P1P1GLSEnum: return NUMNODESP1; - case MINIcondensedEnum: return NUMNODESP1; - case MINIEnum: return NUMNODESP1b; - case TaylorHoodEnum: return NUMNODESP2; - default: _error_("Element type "<element_type)<<" not supported yet"); - } - - return -1; -} -/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 15687) +++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h (revision 15688) @@ -71,5 +71,4 @@ int NumberofNodes(void); int NumberofNodesVelocity(void); - int NumberofNodesVelocityFinal(void); int NumberofNodesPressure(void); }; Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 15688) @@ -833,5 +833,5 @@ /*intermediary: */ - int i; + int i,M,N; int configuration_type; Element *element = NULL; @@ -846,5 +846,5 @@ /*Display message*/ - if(VerboseModule()) _printf0_(" Generating matrices\n"); + if(VerboseModule()) _printf0_(" Generating matrices"); /*retrive parameters: */ @@ -877,4 +877,10 @@ /*Allocate stiffness matrices and load vector*/ this->AllocateSystemMatrices(&Kff,&Kfs,&df,&pf); + + /*Display size*/ + if(VerboseModule()){ + Kff->GetSize(&M,&N); + _printf0_(" (Kff stiffness matrix size: "<nrows==this->ncols && bsize>0 && bsizencols && this->values); + + /*Intermediaries*/ + int counter,i,j,isize; + IssmDouble *Kii = NULL; + IssmDouble *Kib = NULL; + IssmDouble *Kbi = NULL; + IssmDouble *Kbb = NULL; + IssmDouble *Kbbinv = NULL; + IssmDouble *Ktemp = NULL; + int *iindices = NULL; + bool flag; + + /*Get new sizes and indices*/ + isize = this->nrows - bsize; + iindices = xNew(isize); + counter = 0; + for(i=0;inrows;i++){ + flag = true; + for(j=0;j(isize*isize); + Kib = xNew(isize*bsize); + Kbi = xNew(bsize*isize); + Kbb = xNew(bsize*bsize); + for(i=0;ivalues[iindices[i]*this->ncols + iindices[j]]; + for(i=0;ivalues[iindices[i]*this->ncols + bindices[j]]; + for(i=0;ivalues[bindices[i]*this->ncols + iindices[j]]; + for(i=0;ivalues[bindices[i]*this->ncols + bindices[j]]; + + /*Invert Kbb*/ + Kbbinv = xNew(bsize*bsize); + switch(bsize){ + case 1: + Kbbinv[0] = 1./Kbb[0]; + break; + case 2: + Matrix2x2Invert(Kbbinv,Kbb); + break; + case 3: + Matrix3x3Invert(Kbbinv,Kbb); + break; + default: + MatrixInverse(Kbbinv,bsize,bsize,NULL,0,NULL); + break; + } + + /*Calculate Kib inv(Kbb) Kbi*/ + Ktemp = xNew(isize*isize); + TripleMultiply(Kib,isize,bsize,0, Kbbinv,bsize,bsize,0, Kbi,bsize,isize,0, Ktemp,0); + + /*New Ke*/ + for(i=0;inrows*this->ncols;i++) this->values[i]=0.; + for(i=0;ivalues[iindices[i]*this->ncols + iindices[j]] = Ktemp[i*isize+j]; + } + } + + /*Clean up and return*/ + xDelete(Kii); + xDelete(Kib); + xDelete(Kbi); + xDelete(Kbb); + xDelete(Kbbinv); + xDelete(Ktemp); + xDelete(iindices); + return; +} +/*}}}*/ /*FUNCTION ElementMatrix::Echo{{{*/ void ElementMatrix::Echo(void){ Index: /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h (revision 15687) +++ /issm/trunk-jpl/src/c/classes/matrix/ElementMatrix.h (revision 15688) @@ -10,5 +10,4 @@ /*Headers:*/ -/*{{{*/ #include "../../datastructures/datastructures.h" #include "../../toolkits/toolkits.h" @@ -17,5 +16,4 @@ template class Matrix; class Parameters; -/*}}}*/ class ElementMatrix{ @@ -51,5 +49,5 @@ int* col_sglobaldoflist; - /*ElementMatrix constructors, destructors {{{*/ + /*ElementMatrix constructors, destructors*/ ElementMatrix(); ElementMatrix(ElementMatrix* Ke); @@ -58,14 +56,14 @@ ElementMatrix(Node** nodes,int numnodes,Parameters* parameters,int approximation=NoneApproximationEnum); ~ElementMatrix(); - /*}}}*/ - /*ElementMatrix specific routines {{{*/ + + /*ElementMatrix specific routines*/ void AddToGlobal(Matrix* Kff, Matrix* Kfs); void AddToGlobal(Matrix* Jff); void Echo(void); void CheckConsistency(void); + void StaticCondensation(int numindices,int* indices); void Transpose(void); void Init(ElementMatrix* Ke); void SetDiag(IssmDouble scalar); - /*}}}*/ }; #endif //#ifndef _ELEMENT_MATRIX_H_ Index: /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/classes/matrix/ElementVector.cpp (revision 15688) @@ -274,2 +274,105 @@ } /*}}}*/ +/*FUNCTION ElementVector::StaticCondensation{{{*/ +void ElementVector::StaticCondensation(ElementMatrix* Ke,int bsize,int* bindices){ + /* + * | Kii Kib | | Ui | |Fi| + * | Kbi Kbb | | Ub | = |Fb| + * + * Kii Ui + Kib Ub = Fi + * Kbi Ui + Kbb Ub = Fb + * + * We want to remove Ub from the equation: + * + * Kii Ui + Kib inv(Kbb) (Fb - Kbi Ui) = Fi + * + * which gives: + * + * (Kii - Kib inv(Kbb) Kbi) Ui = Fi - Kib inv(Kbb) Fb + */ + + /*Checks in debugging mode*/ + _assert_(bsize>0 && bsizenrows && this->values && Ke); + _assert_(Ke->nrows==Ke->ncols && this->nrows==Ke->nrows); + + /*Intermediaries*/ + int counter,i,j,isize; + IssmDouble *Fb = NULL; + IssmDouble *Fi = NULL; + IssmDouble *Kib = NULL; + IssmDouble *Kbb = NULL; + IssmDouble *Kbbinv = NULL; + IssmDouble *Ftemp = NULL; + int *iindices = NULL; + bool flag; + + /*Get new sizes and indices*/ + isize = this->nrows - bsize; + iindices = xNew(isize); + counter = 0; + for(i=0;inrows;i++){ + flag = true; + for(j=0;j(isize*bsize); + Kbb = xNew(bsize*bsize); + Fb = xNew(bsize); + Fi = xNew(isize); + for(i=0;ivalues[iindices[i]*Ke->ncols + bindices[j]]; + for(i=0;ivalues[bindices[i]*Ke->ncols + bindices[j]]; + for(i=0;ivalues[bindices[i]]; + for(i=0;ivalues[iindices[i]]; + + /*Invert Kbb*/ + Kbbinv = xNew(bsize*bsize); + switch(bsize){ + case 1: + Kbbinv[0] = 1./Kbb[0]; + break; + case 2: + Matrix2x2Invert(Kbbinv,Kbb); + break; + case 3: + Matrix3x3Invert(Kbbinv,Kbb); + break; + default: + MatrixInverse(Kbbinv,bsize,bsize,NULL,0,NULL); + break; + } + + /*Calculate Kib inv(Kbb) Fb*/ + Ftemp = xNew(isize); + TripleMultiply(Kib,isize,bsize,0, Kbbinv,bsize,bsize,0, Fb,bsize,1,0, Ftemp,0); + + /*New Pe*/ + for(i=0;inrows;i++) this->values[i]=0.; + for(i=0;ivalues[iindices[i]] = Ftemp[i]; + } + + /*Clean up and return*/ + xDelete(Kib); + xDelete(Kbb); + xDelete(Kbbinv); + xDelete(Fb); + xDelete(Fi); + xDelete(Ftemp); + xDelete(iindices); + return; +} +/*}}}*/ Index: /issm/trunk-jpl/src/c/classes/matrix/ElementVector.h =================================================================== --- /issm/trunk-jpl/src/c/classes/matrix/ElementVector.h (revision 15687) +++ /issm/trunk-jpl/src/c/classes/matrix/ElementVector.h (revision 15688) @@ -10,5 +10,4 @@ /*Headers:*/ -/*{{{*/ #include "../../datastructures/datastructures.h" #include "../../toolkits/toolkits.h" @@ -17,22 +16,21 @@ template class Vector; class Parameters; -/*}}}*/ +class ElementMatrix; class ElementVector{ public: - - int nrows; - IssmDouble* values; + int nrows; + IssmDouble* values; //gset - int* gglobaldoflist; + int* gglobaldoflist; //fset - int fsize; - int* flocaldoflist; - int* fglobaldoflist; + int fsize; + int* flocaldoflist; + int* fglobaldoflist; - /*ElementVector constructors, destructors {{{*/ + /*ElementVector constructors, destructors*/ ElementVector(); ElementVector(ElementVector* pe1,ElementVector* pe2); @@ -40,6 +38,6 @@ ElementVector(Node** nodes,int numnodes,Parameters* parameters,int approximation=NoneApproximationEnum); ~ElementVector(); - /*}}}*/ - /*ElementVector specific routines {{{*/ + + /*ElementVector specific routines*/ void AddToGlobal(Vector* pf); void InsertIntoGlobal(Vector* pf); @@ -48,5 +46,5 @@ void Init(ElementVector* pe); void SetValue(IssmDouble scalar); - /*}}}*/ + void StaticCondensation(ElementMatrix* Ke,int numindices,int* indices); }; #endif //#ifndef _ELEMENT_VECTOR_H_ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 15688) @@ -130,14 +130,21 @@ case MINIcondensedEnum: _assert_(approximation==FSApproximationEnum); - /*P1 velocity (bubble statically condensed)*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum)); - } - } - /*P1 pressure*/ - for(i=0;inumberofvertices;i++){ - if(iomodel->my_vertices[i]){ - nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,i,iomodel,analysis,FSpressureEnum)); + /*P1+ velocity (bubble statically condensed)*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i,iomodel,analysis,FSvelocityEnum)); + } + } + for(i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,analysis,FSvelocityEnum); + node->Deactivate(); + nodes->AddObject(node); + } + } + /*P1 pressure*/ + for(i=0;inumberofvertices;i++){ + if(iomodel->my_vertices[i]){ + nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,i,iomodel,analysis,FSpressureEnum)); } } Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp (revision 15688) @@ -98,5 +98,5 @@ case 0 : finiteelement = P1Enum; break;//P1P1 case 1 : finiteelement = P1Enum; break;//P1P1GSL - case 2 : finiteelement = P1Enum; break;//MINIcondensed + case 2 : finiteelement = P1bubbleEnum; break;//MINIcondensed case 3 : finiteelement = P1bubbleEnum; break;//MINI case 4 : finiteelement = P2Enum; break;//TaylorHood (P2P1) Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp (revision 15688) @@ -70,5 +70,4 @@ } } - iomodel->FetchData(9,MeshVertexonbedEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum, MaskVertexongroundediceEnum,MaskVertexonfloatingiceEnum,MaskVertexonwaterEnum,FlowequationVertexEquationEnum,DiagnosticReferentialEnum); Index: /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 15687) +++ /issm/trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp (revision 15688) @@ -21,7 +21,7 @@ int idima,idimb,idimc,idimd; IssmDouble *dtemp; + _assert_(a && b && c && d); /* set up dimensions for triple product */ - if (!itrna){ idima=nrowa; @@ -83,4 +83,6 @@ int i,j,k,ipta,iptb,iptc; int nrowc,ncolc,iinca,jinca,iincb,jincb,ntrma,ntrmb,nterm; + + _assert_(a && b && c); /* set up dimensions and increments for matrix a */ @@ -326,5 +328,5 @@ /*Compute determinant*/ Matrix2x2Determinant(&det,A); - if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller that machine epsilon"); + if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon"); /*Compute invert*/ @@ -351,5 +353,5 @@ /*Compute determinant*/ Matrix3x3Determinant(&det,A); - if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller that machine epsilon"); + if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon"); /*Compute invert*/