[cig-commits] r6490 - in short/3D/PyLith/trunk: . libsrc/meshio
modulesrc/feassemble modulesrc/topology pylith/feassemble
pylith/problems unittests/pytests/meshio/data
brad at geodynamics.org
brad at geodynamics.org
Sat Mar 31 14:50:13 PDT 2007
Author: brad
Date: 2007-03-31 14:50:12 -0700 (Sat, 31 Mar 2007)
New Revision: 6490
Modified:
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src
short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt
Log:
Added methods for creating PETSc matrix.
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/TODO 2007-03-31 21:50:12 UTC (rev 6490)
@@ -7,19 +7,17 @@
add check to material::initialize material dimension must match cell dimension
-0. Update Material tests. Add unit tests for ElasticIsotropic1D,
-ElasticPlaneStrain, ElasticPlaneStress.
+0. Update unit testing.
+ a. Update Material tests.
+ b. Add unit tests for ElasticIsotropic1D, ElasticPlaneStrain,
+ ElasticPlaneStress.
+ c. Update MeshIOAscii tests (groups of vertices and cells)
+
1. Finish implementing ExplicitElasticity
- a. C++
- b. Python object
- c. bindings
d. unit tests at C++ level
e. unit tests at Python level
- Status: Started on a.
- Only implement nonlumped and let PETSc "solver" do the lumping.
-
2. Finish implementing Python Formulation
(e.g., Explicit) with initialization of solid element families.
a. Python object
@@ -29,9 +27,6 @@
SECONDARY PRIORITIES
======================================================================
-0. Adjust MeshIOAscii format
- b. Add group of vertices
-
1. Implement MeshIOHDF5 & HDF5
a. C++ objects
b. unit tests at C++ level
Modified: short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/libsrc/meshio/MeshIOAscii.cc 2007-03-31 21:50:12 UTC (rev 6490)
@@ -20,7 +20,8 @@
#include <assert.h> // USES assert()
#include <iomanip> // USES setw(), setiosflags(), resetiosflags()
-const char *pylith::meshio::MeshIOAscii::groupTypeNames[] = {"vertices", "cells"};
+const char* pylith::meshio::MeshIOAscii::groupTypeNames[] =
+ {"vertices", "cells"};
// ----------------------------------------------------------------------
// Constructor
@@ -155,9 +156,6 @@
_writeGroup(fileout, name->c_str());
}
- // LOOP OVER GROUPS
- // _writeGroup(fileout, mesh, nameIter->c_str());
-
fileout << "}\n";
fileout.close();
} // write
@@ -401,6 +399,8 @@
int& numPoints,
int *points[]) const
{ // _readGroup
+ assert(0 != points);
+
int* indices = 0; // Indices of entities in group
std::string token;
@@ -435,9 +435,8 @@
} // if
filein.ignore(maxIgnore, '{');
- delete[] indices; indices = new int[numPoints];
- assert(0 != indices);
- for (int i = 0; i < numPoints; ++i)
+ delete[] indices; indices = (numPoints > 0) ? new int[numPoints] : 0;
+ for (int i=0; i < numPoints; ++i)
filein >> indices[i];
filein.ignore(maxIgnore, '}');
} else {
@@ -458,7 +457,7 @@
pylith::meshio::MeshIOAscii::_writeGroup(std::ostream& fileout,
const char* name) const
{ // _writeGroup
- pylith::meshio::MeshIO::PointType type;
+ pylith::meshio::MeshIO::PointType type;
int numPoints;
int *points;
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2007-03-31 21:50:12 UTC (rev 6490)
@@ -23,6 +23,8 @@
#include "pylith/feassemble/IntegratorExplicit.hh"
#include "pylith/feassemble/ExplicitElasticity.hh"
+#include "pylith/utils/petscfwd.h"
+
#include <assert.h>
#include <stdexcept>
#include <Python.h>
@@ -495,7 +497,40 @@
"""
Compute matrix (A) associated with operator.
"""
- print "WARNING: IntegratorExplicit::integrateJacobian not implemented."
+ # create shim for method 'integrateJacobian'
+ #embed{ void IntegratorExplicit_integrateJacobian(void* objVptr, void* matVptr, void* fieldInVptr, void* meshVptr)
+ typedef ALE::Field::Mesh Mesh;
+ typedef ALE::Field::Mesh::real_section_type real_section_type;
+
+ try {
+ assert(0 != objVptr);
+ assert(0 != matVptr);
+ assert(0 != fieldInVptr);
+ assert(0 != meshVptr);
+ ALE::Obj<Mesh>* mesh =
+ (ALE::Obj<Mesh>*) meshVptr;
+ PetscMat* mat = (PetscMat*) matVptr;
+ ALE::Obj<real_section_type>* fieldIn =
+ (ALE::Obj<real_section_type>*) fieldInVptr;
+ ((pylith::feassemble::IntegratorExplicit*) objVptr)->integrateJacobian(
+ mat, *fieldIn, *mesh);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ cdef void* matVptr
+ cdef void* fieldInVptr
+ matVptr = PyCObject_AsVoidPtr(mat)
+ fieldInVptr = PyCObject_AsVoidPtr(fieldIn)
+ IntegratorExplicit_integrateJacobian(self.thisptr, matVptr, fieldInVptr,
+ ptrFromHandle(mesh))
return
Modified: short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/modulesrc/topology/topology.pyxe.src 2007-03-31 21:50:12 UTC (rev 6490)
@@ -12,6 +12,8 @@
#header{
#include <Mesh.hh>
+#include <petscmesh.h>
+#include "pylith/utils/petscfwd.h"
#include <stdexcept>
#include <Python.h>
@@ -43,6 +45,20 @@
MeshPtr_destructor_cpp(obj)
return
+cdef void PetscMat_destructor(void* obj):
+ """
+ Destroy PetscMat.
+ """
+ #embed{ void PetscMat_destructor_cpp(void* objVptr)
+ assert(0 != objVptr);
+ PetscMat* mat = (PetscMat*) objVptr;
+ PetscErrorCode err = MatDestroy(*mat);
+ if (err)
+ throw std::runtime_error("Could not destroy PETSc matrix.");
+ #}embed
+ PetscMat_destructor_cpp(obj)
+ return
+
# ----------------------------------------------------------------------
cdef class Mesh:
@@ -122,6 +138,7 @@
assert(!section.isNull());
section->setFiberDimension((*mesh)->depthStratum(0), fiberDim);
(*mesh)->allocate(section);
+ section->zero();
result = (void*) §ion;
} catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError,
@@ -137,6 +154,43 @@
return PyCObject_FromVoidPtr(ptr, NULL)
+ def createMatrix(self, field):
+ """
+ Create matrix compatible with field.
+ """
+ # create shim for MeshCreateMatrix
+ #embed{ void* Mesh_createMatrix(void* objVptr, void* fieldVptr)
+ typedef ALE::Field::Mesh::real_section_type real_section_type;
+
+ void* result = 0;
+ try {
+ ALE::Obj<ALE::Field::Mesh>* mesh = (ALE::Obj<ALE::Field::Mesh>*) objVptr;
+ assert(0 != mesh);
+ assert(!mesh->isNull());
+ const ALE::Obj<real_section_type>* field =
+ (ALE::Obj<real_section_type>*) fieldVptr;
+ assert(!field->isNull());
+
+ PetscMat* mat = new PetscMat;
+ assert(0 != mat);
+ MeshCreateMatrix(*mesh, *field, MATMPIBAIJ, mat);
+ result = (void*) mat;
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ return result;
+ #}embed
+ cdef void* ptr
+ cdef void* fieldVptr
+ fieldVptr = PyCObject_AsVoidPtr(field)
+ ptr = Mesh_createMatrix(self.thisptr, fieldVptr)
+ return PyCObject_FromVoidPtr(ptr, PetscMat_destructor)
+
+
def _createHandle(self):
"""
Wrap pointer to C++ object in PyCObject.
Modified: short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/pylith/feassemble/IntegratorExplicit.py 2007-03-31 21:50:12 UTC (rev 6490)
@@ -62,7 +62,7 @@
return
- def integrateJacobian(self, jacobian, fieldInT, coords):
+ def integrateJacobian(self, jacobian, fieldInT):
"""
Integrate Jacobian term for dynamic terms for finite-elements.
"""
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2007-03-31 21:50:12 UTC (rev 6490)
@@ -87,15 +87,15 @@
self.integrators.append(integrator)
self._info.log("Creating fields and matrices.")
- # self.jacobian = mesh.cppHandle.getPetscMat()
self.dispT = mesh.cppHandle.createRealSection("dispT", spaceDim)
self.dispTmdt = mesh.cppHandle.createRealSection("dispTmdt", spaceDim)
self.dispTpdt = mesh.cppHandle.createRealSection("dispTpdt", spaceDim)
self.constant = mesh.cppHandle.createRealSection("constant", spaceDim)
+ self.jacobian = mesh.cppHandle.createMatrix(self.constant)
self._info.log("Integrating Jacobian of operator.")
- #for integrator in integrators:
- # integrator.integrateJacobian(self.jacobian, self.dispT)
+ for integrator in self.integrators:
+ integrator.integrateJacobian(self.jacobian, self.dispT)
return
Modified: short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt 2007-03-31 19:54:13 UTC (rev 6489)
+++ short/3D/PyLith/trunk/unittests/pytests/meshio/data/mesh2Din3D.txt 2007-03-31 21:50:12 UTC (rev 6490)
@@ -30,4 +30,20 @@
2 0
}
}
+ group = {
+ type = vertices
+ name = group A
+ count = 4
+ indices = {
+ 1 3 5 7
+ }
+ }
+ group = {
+ type = cells
+ name = group B
+ count = 1
+ indices = {
+ 1
+ }
+ }
}
More information about the cig-commits
mailing list