[cig-commits] r4018 - in
short/3D/PyLith/branches/pylith-0.8/pylith3d: libpylith3d
module pylith3d
knepley at geodynamics.org
knepley at geodynamics.org
Fri Jul 14 10:20:44 PDT 2006
Author: knepley
Date: 2006-07-14 10:20:43 -0700 (Fri, 14 Jul 2006)
New Revision: 4018
Modified:
short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/makemsr.F
short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc
short/3D/PyLith/branches/pylith-0.8/pylith3d/pylith3d/Pylith3d_run.py
Log:
Updates for new Sieve stuff
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/makemsr.F
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/makemsr.F 2006-07-13 23:41:51 UTC (rev 4017)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/makemsr.F 2006-07-14 17:20:43 UTC (rev 4018)
@@ -71,7 +71,7 @@
nmin=10000
nmax=izero
wavg=dble(neq)
- call MatSetFromOptions(A, ierr)
+c call MatSetFromOptions(A, ierr)
call IsCreateStride(MPI_COMM_SELF, neq, 1, 1, rowSizes, ierr)
call ISGetIndices(rowSizes, row_array, row_offset, ierr)
do i=1,neq
@@ -89,9 +89,9 @@
nmin=nmin+ione
nmax=nmax+ione
if(neq.ne.izero) wavg=wavg/dble(neq)
-cFIX call MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, ierr)
- call MatSeqAIJSetPreallocation(A, 0, row_array(row_offset+1),
- & ierr)
+c call MatSetOption(A, MAT_NEW_NONZERO_ALLOCATION_ERR, ierr)
+c call MatSeqAIJSetPreallocation(A, 0, row_array(row_offset+1),
+c & ierr)
call ISRestoreIndices(rowSizes, row_array, row_offset, ierr)
call ISDestroy(rowSizes, ierr)
return
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc 2006-07-13 23:41:51 UTC (rev 4017)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc 2006-07-14 17:20:43 UTC (rev 4018)
@@ -313,7 +313,7 @@
// return
PyObject *pyMesh = PyCObject_FromVoidPtr(mesh.ptr(), NULL);
- mesh.int_allocator.del(mesh.refCnt);
+ mesh.int_allocator->del(mesh.refCnt);
mesh.refCnt = NULL;
return Py_BuildValue("sN", meshOutputFile, pyMesh);
}
@@ -384,6 +384,9 @@
ierr = PetscObjectContainerDestroy(c);
ierr = PetscObjectCompose((PetscObject) sol, "injection", (PetscObject) injection);
+ ierr = MatSetFromOptions(A);
+ ierr = preallocateMatrix(A, mesh, mesh->getField("displacement"));
+
journal::debug_t debug("pylith3d");
debug
<< journal::at(__HERE__)
@@ -546,7 +549,7 @@
PetscViewerPopFormat(viewer);
PetscViewerDestroy(viewer);
- m.int_allocator.del(m.refCnt);
+ m.int_allocator->del(m.refCnt);
m.refCnt = NULL;
journal::debug_t debug("pylith3d");
@@ -560,6 +563,112 @@
return Py_None;
}
+char pypylith3d_interpolatePoints__doc__[] = "";
+char pypylith3d_interpolatePoints__name__[] = "interpolatePoints";
+
+PyObject * pypylith3d_interpolatePoints(PyObject *, PyObject *args)
+{
+ PyObject *pyMesh, *pySol, *pyPoints,*pyValues;
+
+ int ok = PyArg_ParseTuple(args, "OOO:outputMesh", &pyMesh, &pySol, &pyPoints);
+ if (!ok) {
+ return 0;
+ }
+
+ ALE::Mesh *mesh = (ALE::Mesh *) PyCObject_AsVoidPtr(pyMesh);
+ Vec sol = (Vec) PyCObject_AsVoidPtr(pySol);
+ double *points;
+
+ // Convert Numeric matrix to C array
+
+ ALE::Obj<ALE::Mesh> m(mesh);
+
+ // Injection Vec in to Field
+ ALE::Obj<ALE::Mesh::field_type> displacement = m->getField("displacement");
+ ALE::Mesh::field_type::patch_type patch;
+ Vec l;
+ VecScatter injection;
+
+ VecCreateSeqWithArray(PETSC_COMM_SELF, displacement->getSize(patch), displacement->restrict(patch), &l);
+ PetscObjectQuery((PetscObject) sol, "injection", (PetscObject *) &injection);
+ VecScatterBegin(sol, l, INSERT_VALUES, SCATTER_REVERSE, injection);
+ VecScatterEnd(sol, l, INSERT_VALUES, SCATTER_REVERSE, injection);
+ VecDestroy(l);
+
+ // Create complete field by adding BC
+ ALE::Obj<ALE::Mesh::field_type> full_displacement = m->getField("full_displacement");
+ ALE::Obj<ALE::Mesh::field_type::order_type::baseSequence> patches = displacement->getPatches();
+ ALE::Obj<ALE::Mesh::foliation_type> boundaries = m->getBoundaries();
+
+ // This is wrong if the domain changes
+ if (!full_displacement->getGlobalOrder()) {
+ for(ALE::Mesh::field_type::order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
+ full_displacement->setPatch(displacement->getPatch(*p_iter), *p_iter);
+ full_displacement->setFiberDimensionByDepth(*p_iter, 0, 3);
+ }
+ full_displacement->orderPatches();
+ full_displacement->createGlobalOrder();
+ }
+ for(ALE::Mesh::field_type::order_type::baseSequence::iterator p_iter = patches->begin(); p_iter != patches->end(); ++p_iter) {
+ ALE::Obj<ALE::Mesh::field_type::order_type::coneSequence> elements = full_displacement->getPatch(*p_iter);
+
+ for(ALE::Mesh::field_type::order_type::coneSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
+ const int dim = displacement->getIndex(*p_iter, *e_iter).index;
+ const double *array = displacement->restrict(*p_iter, *e_iter);
+ int v = 0;
+ double values[3];
+
+ for(int c = 0; c < 3; c++) {
+ ALE::Mesh::foliation_type::patch_type bPatch(*p_iter, c+1);
+ const ALE::Mesh::foliation_type::index_type& idx = boundaries->getIndex(bPatch, *e_iter);
+
+ if (idx.index > 0) {
+ values[c] = 0.0;
+ } else if (dim > 0) {
+ values[c] = array[v++];
+ }
+ }
+ if (v != dim) {
+ std::cout << "ERROR: Invalid size " << v << " used for " << *e_iter << " with index " << displacement->getIndex(*p_iter, *e_iter) << std::endl;
+ }
+ full_displacement->updateAdd(*p_iter, *e_iter, values);
+ }
+ }
+ ALE::Obj<ALE::Mesh::foliation_type::order_type::baseSequence> bdPatches = boundaries->getPatches();
+
+ for(ALE::Mesh::foliation_type::order_type::baseSequence::iterator p_iter = bdPatches->begin(); p_iter != bdPatches->end(); ++p_iter) {
+ ALE::Obj<ALE::Mesh::foliation_type::order_type::coneSequence> elements = boundaries->getPatch(*p_iter);
+
+ for(ALE::Mesh::foliation_type::order_type::coneSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
+ const double *array = full_displacement->restrict((*p_iter).first, *e_iter);
+ double values[3];
+
+ if (boundaries->getIndex(*p_iter, *e_iter).index > 0) {
+ for(int c = 0; c < 3; c++) {
+ values[c] = array[c];
+ }
+ values[(*p_iter).second-1] = boundaries->restrict(*p_iter, *e_iter)[0];
+ full_displacement->update((*p_iter).first, *e_iter, values);
+ }
+ }
+ }
+
+ //interpolatePoints(m, full_displacement, points, values);
+
+ // Convert C array to Numeric matrix
+
+ m.int_allocator->del(m.refCnt);
+ m.refCnt = NULL;
+
+ journal::debug_t debug("pylith3d");
+ debug
+ << journal::at(__HERE__)
+ << "Interpolated points"
+ << journal::endl;
+
+ return Py_BuildValue("N", pyValues);
+}
+
// Scan boundary conditions
char pypylith3d_scan_bc__doc__[] = "";
Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/pylith3d/Pylith3d_run.py
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/pylith3d/Pylith3d_run.py 2006-07-13 23:41:51 UTC (rev 4017)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/pylith3d/Pylith3d_run.py 2006-07-14 17:20:43 UTC (rev 4018)
@@ -178,6 +178,103 @@
return
+ def solveElastic(self):
+ import pylith3d
+ pylith3d.elastc(
+ self.A,self.rhs,self.sol, # sparse
+ self.pointerToBextern, # force
+ self.pointerToBtraction,
+ self.pointerToBgravity,
+ self.pointerToBconcForce,
+ self.pointerToBintern,
+ self.pointerToBresid,
+ self.pointerToBwink,
+ self.pointerToBwinkx,
+ self.pointerToDispVec,
+ self.pointerToDprev,
+ self.pointerToListArrayNforce,
+ self.pointerToListArrayGrav,
+ self.pointerToX, # global
+ self.pointerToD,
+ self.pointerToDeld,
+ self.pointerToDcur,
+ self.pointerToId,
+ self.pointerToIwink,
+ self.pointerToWink,
+ self.pointerToListArrayNsysdat,
+ self.pointerToListArrayIddmat,
+ self.pointerToIbond, # BC
+ self.pointerToBond,
+ self.pointerToDx, # slip
+ self.pointerToDeldx,
+ self.pointerToDxcur,
+ self.pointerToDiforc,
+ self.pointerToIdx,
+ self.pointerToIwinkx,
+ self.pointerToWinkx,
+ self.pointerToIdslp,
+ self.pointerToIpslp,
+ self.pointerToIdhist,
+ self.pointerToFault, # fault
+ self.pointerToNfault,
+ self.pointerToDfault,
+ self.pointerToTfault,
+ self.pointerToS, # stiff
+ self.pointerToStemp,
+ self.pointerToState, # element
+ self.pointerToDstate,
+ self.pointerToState0,
+ self.pointerToDmat,
+ self.pointerToIens,
+ self.pointerToLm,
+ self.pointerToLmx,
+ self.pointerToLmf,
+ self.pointerToIvfamily,
+ self.pointerToListArrayNpar,
+ self.pointerToIelindx,
+ self.pointerToIelno, # traction
+ self.pointerToIside,
+ self.pointerToIhistry,
+ self.pointerToPres,
+ self.pointerToPdir,
+ self.pointerToListArrayPropertyList, # material
+ self.pointerToMaterialModelInfo,
+ self.pointerToGauss, # eltype
+ self.pointerToSh,
+ self.pointerToShj,
+ self.pointerToListArrayElementTypeInfo,
+ self.pointerToHistry, # timdat
+ self.pointerToListArrayRtimdat,
+ self.pointerToListArrayNtimdat,
+ self.pointerToListArrayNvisdat,
+ self.pointerToMaxstp,
+ self.pointerToDelt,
+ self.pointerToAlfa,
+ self.pointerToMaxit,
+ self.pointerToNtdinit,
+ self.pointerToLgdef,
+ self.pointerToUtol,
+ self.pointerToFtol,
+ self.pointerToEtol,
+ self.pointerToItmax,
+ self.pointerToListArrayRgiter, # stresscmp
+ self.pointerToSkew, # skew
+ self.pointerToListArrayNcodat, # ioinfo
+ self.pointerToListArrayNunits,
+ self.pointerToListArrayNprint,
+ self.pointerToIstatout,
+ self.pointerToNstatout,
+ self.asciiOutputFile, # files
+ self.plotOutputFile,
+ self.ucdOutputRoot,
+ self.elasticStage, # PETSc logging
+ self.iterateEvent)
+ return
+
+ def interpolatePoints(self, points):
+ import pylith3d
+ return pylith3d.(self.mesh, self.sol, points)
+
def run(self):
import pylith3d
@@ -292,97 +389,7 @@
self.iterateEvent)
# Perform elastic solution, if requested.
-
- pylith3d.elastc(
- self.A,self.rhs,self.sol, # sparse
- self.pointerToBextern, # force
- self.pointerToBtraction,
- self.pointerToBgravity,
- self.pointerToBconcForce,
- self.pointerToBintern,
- self.pointerToBresid,
- self.pointerToBwink,
- self.pointerToBwinkx,
- self.pointerToDispVec,
- self.pointerToDprev,
- self.pointerToListArrayNforce,
- self.pointerToListArrayGrav,
- self.pointerToX, # global
- self.pointerToD,
- self.pointerToDeld,
- self.pointerToDcur,
- self.pointerToId,
- self.pointerToIwink,
- self.pointerToWink,
- self.pointerToListArrayNsysdat,
- self.pointerToListArrayIddmat,
- self.pointerToIbond, # BC
- self.pointerToBond,
- self.pointerToDx, # slip
- self.pointerToDeldx,
- self.pointerToDxcur,
- self.pointerToDiforc,
- self.pointerToIdx,
- self.pointerToIwinkx,
- self.pointerToWinkx,
- self.pointerToIdslp,
- self.pointerToIpslp,
- self.pointerToIdhist,
- self.pointerToFault, # fault
- self.pointerToNfault,
- self.pointerToDfault,
- self.pointerToTfault,
- self.pointerToS, # stiff
- self.pointerToStemp,
- self.pointerToState, # element
- self.pointerToDstate,
- self.pointerToState0,
- self.pointerToDmat,
- self.pointerToIens,
- self.pointerToLm,
- self.pointerToLmx,
- self.pointerToLmf,
- self.pointerToIvfamily,
- self.pointerToListArrayNpar,
- self.pointerToIelindx,
- self.pointerToIelno, # traction
- self.pointerToIside,
- self.pointerToIhistry,
- self.pointerToPres,
- self.pointerToPdir,
- self.pointerToListArrayPropertyList, # material
- self.pointerToMaterialModelInfo,
- self.pointerToGauss, # eltype
- self.pointerToSh,
- self.pointerToShj,
- self.pointerToListArrayElementTypeInfo,
- self.pointerToHistry, # timdat
- self.pointerToListArrayRtimdat,
- self.pointerToListArrayNtimdat,
- self.pointerToListArrayNvisdat,
- self.pointerToMaxstp,
- self.pointerToDelt,
- self.pointerToAlfa,
- self.pointerToMaxit,
- self.pointerToNtdinit,
- self.pointerToLgdef,
- self.pointerToUtol,
- self.pointerToFtol,
- self.pointerToEtol,
- self.pointerToItmax,
- self.pointerToListArrayRgiter, # stresscmp
- self.pointerToSkew, # skew
- self.pointerToListArrayNcodat, # ioinfo
- self.pointerToListArrayNunits,
- self.pointerToListArrayNprint,
- self.pointerToIstatout,
- self.pointerToNstatout,
- self.asciiOutputFile, # files
- self.plotOutputFile,
- self.ucdOutputRoot,
- self.elasticStage, # PETSc logging
- self.iterateEvent)
-
+ self.solveElastic()
pylith3d.outputMesh(self.fileRoot, self.mesh, self.sol)
# Perform time-dependent solution, if requested.
More information about the cig-commits
mailing list