[cig-commits] r20045 - short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Sat May 5 10:49:18 PDT 2012
Author: brad
Date: 2012-05-05 10:49:17 -0700 (Sat, 05 May 2012)
New Revision: 20045
Modified:
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.hh
Log:
Add setting null space in Solver for field split (from Matt).
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.cc 2012-05-05 00:13:07 UTC (rev 20044)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.cc 2012-05-05 17:49:17 UTC (rev 20045)
@@ -168,20 +168,87 @@
bool separateComponents = formulation->splitFieldComponents();
if (separateComponents) {
constructFieldSplit(solutionSection, PETSC_DETERMINE, PETSC_NULL, PETSC_NULL,
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
- solutionSection),
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionSection), PETSC_NULL,
solution.vector(), *pc);
} else {
int numFields[2] = {spaceDim, (numSpaces > spaceDim) ? 1 : 0};
+ MatNullSpace nullsp[2] = {PETSC_NULL, PETSC_NULL};
int* fields = new int[numSpaces];
- for(PetscInt f=0; f < numSpaces; ++f) {
+
+ /* Create rigid body null space */
+ const ALE::Obj<RealSection>& coordinatesSection = sieveMesh->getRealSection("coordinates");
+ assert(!coordinatesSection.isNull());
+ const ALE::Obj<SieveMesh::label_sequence>& vertices = sieveMesh->depthStratum(0);
+ assert(!vertices.isNull());
+ const SieveMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveMesh::label_sequence::iterator verticesEnd = vertices->end();
+ PetscInt dim = spaceDim;
+ Vec mode[6];
+
+ if (dim > 1) {
+ PetscInt n;
+ err = VecGetLocalSize(solution.vector(), &n);CHECK_PETSC_ERROR(err);
+ const int m = (dim * (dim + 1)) / 2;
+ err = VecCreate(sieveMesh->comm(), &mode[0]);CHECK_PETSC_ERROR(err);
+ err = VecSetSizes(mode[0], n, PETSC_DETERMINE);CHECK_PETSC_ERROR(err);
+ err = VecSetUp(mode[0]);CHECK_PETSC_ERROR(err);
+ for (int i = 1; i < m; ++i) {
+ err = VecDuplicate(mode[0], &mode[i]);CHECK_PETSC_ERROR(err);
+ } // for
+ // :KLUDGE: Assume P1
+ for (int d=0; d < dim; ++d) {
+ PetscScalar values[3] = {0.0, 0.0, 0.0};
+
+ values[d] = 1.0;
+ solutionSection->zero();
+ for(SieveMesh::label_sequence::iterator v_iter=verticesBegin; v_iter != verticesEnd; ++v_iter) {
+ solutionSection->updatePoint(*v_iter, values);
+ } // for
+ solution.scatterSectionToVector();
+ err = VecCopy(solution.vector(), mode[d]);CHECK_PETSC_ERROR(err);
+ } // for
+ for (int d = dim; d < m; ++d) {
+ PetscInt k = (dim > 2) ? d - dim : d;
+
+ solutionSection->zero();
+ for (SieveMesh::label_sequence::iterator v_iter=verticesBegin; v_iter != verticesEnd; ++v_iter) {
+ PetscScalar values[3] = {0.0, 0.0, 0.0};
+ const PylithScalar* coords = coordinatesSection->restrictPoint(*v_iter);
+
+ for (int i=0; i < dim; ++i) {
+ for (int j=0; j < dim; ++j) {
+ values[j] += _epsilon(i, j, k)*coords[i];
+ } // for
+ } // for
+ solutionSection->updatePoint(*v_iter, values);
+ }
+ solution.scatterSectionToVector();
+ err = VecCopy(solution.vector(), mode[d]);CHECK_PETSC_ERROR(err);
+ } // for
+ for (int i=0; i < dim; ++i) {
+ err = VecNormalize(mode[i], PETSC_NULL);CHECK_PETSC_ERROR(err);
+ } // for
+ /* Orthonormalize system */
+ for (int i = dim; i < m; ++i) {
+ PetscScalar dots[6];
+
+ err = VecMDot(mode[i], i, mode, dots);CHECK_PETSC_ERROR(err);
+ for(int j=0; j < i; ++j) dots[j] *= -1.0;
+ err = VecMAXPY(mode[i], i, dots, mode);CHECK_PETSC_ERROR(err);
+ err = VecNormalize(mode[i], PETSC_NULL);CHECK_PETSC_ERROR(err);
+ } // for
+ err = MatNullSpaceCreate(sieveMesh->comm(), PETSC_FALSE, m, mode, &nullsp[0]);CHECK_PETSC_ERROR(err);
+ for(int i=0; i< m; ++i) {err = VecDestroy(&mode[i]);CHECK_PETSC_ERROR(err);}
+ } // if
+
+ for (PetscInt f=0; f < numSpaces; ++f) {
fields[f] = f;
} // for
constructFieldSplit(solutionSection, 2, numFields, fields,
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh,
-"default", solutionSection),
- solution.vector(), *pc);
- delete[] fields;
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", solutionSection), nullsp,
+ solution.vector(), *pc);
+ err = MatNullSpaceDestroy(&nullsp[0]);CHECK_PETSC_ERROR(err);
+ delete [] fields;
} // if/else
if (formulation->splitFields() &&
@@ -237,5 +304,64 @@
} // if
} // _setupFieldSplit
+// ----------------------------------------------------------------------
+int
+pylith::problems::Solver::_epsilon(int i,
+ int j,
+ int k)
+{ // _epsilon
+ switch(i) {
+ case 0:
+ switch(j) {
+ case 0: return 0;
+ case 1:
+ switch(k) {
+ case 0: return 0;
+ case 1: return 0;
+ case 2: return 1;
+ }
+ case 2:
+ switch(k) {
+ case 0: return 0;
+ case 1: return -1;
+ case 2: return 0;
+ }
+ }
+ case 1:
+ switch(j) {
+ case 0:
+ switch(k) {
+ case 0: return 0;
+ case 1: return 0;
+ case 2: return -1;
+ }
+ case 1: return 0;
+ case 2:
+ switch(k) {
+ case 0: return 1;
+ case 1: return 0;
+ case 2: return 0;
+ }
+ }
+ case 2:
+ switch(j) {
+ case 0:
+ switch(k) {
+ case 0: return 0;
+ case 1: return 1;
+ case 2: return 0;
+ }
+ case 1:
+ switch(k) {
+ case 0: return -1;
+ case 1: return 0;
+ case 2: return 0;
+ }
+ case 2: return 0;
+ }
+ }
+ return 0;
+} // _epsilon
+
// End of file
Modified: short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.hh
===================================================================
--- short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.hh 2012-05-05 00:13:07 UTC (rev 20044)
+++ short/3D/PyLith/branches/v1.7-trunk/libsrc/pylith/problems/Solver.hh 2012-05-05 17:49:17 UTC (rev 20045)
@@ -91,6 +91,14 @@
const topology::Jacobian& jacobian,
const topology::SolutionFields& fields);
+ /**
+ */
+ static
+ int
+ _epsilon(int i,
+ int j,
+ int k);
+
// PROTECTED MEMBERS ////////////////////////////////////////////////////
protected :
More information about the CIG-COMMITS
mailing list