[cig-commits] r4554 - short/3D/PyLith/trunk/playpen/integrate/src
knepley at geodynamics.org
knepley at geodynamics.org
Fri Sep 15 10:38:57 PDT 2006
Author: knepley
Date: 2006-09-15 10:38:57 -0700 (Fri, 15 Sep 2006)
New Revision: 4554
Modified:
short/3D/PyLith/trunk/playpen/integrate/src/Makefile.am
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc
Log:
Now we have a laplacian application routine
Modified: short/3D/PyLith/trunk/playpen/integrate/src/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/Makefile.am 2006-09-15 16:03:37 UTC (rev 4553)
+++ short/3D/PyLith/trunk/playpen/integrate/src/Makefile.am 2006-09-15 17:38:57 UTC (rev 4554)
@@ -9,8 +9,8 @@
#
# ----------------------------------------------------------------------
#
-MESHIO_INCLUDE = -I../../meshio/src
-MESHIO_LIB = ../../meshio/src/MeshIO.o ../../meshio/src/MeshIOAscii.o
+CIG_INCLUDE = -I/PETSc3/geoframe/packages/pylith-new/include
+CIG_LIB = -L/PETSc3/geoframe/packages/pylith-new/lib -lpylithmeshio
bin_PROGRAMS = testintegrate
@@ -21,9 +21,9 @@
noinst_HEADERS = \
elemVector.hh
-testintegrate_LDADD = $(MESHIO_LIB) $(PETSC_LIB)
+testintegrate_LDADD = $(CIG_LIB) $(PETSC_LIB)
-INCLUDES = $(PETSC_INCLUDE) $(MESHIO_INCLUDE)
+INCLUDES = $(CIG_INCLUDE) $(PETSC_INCLUDE)
# version
# $Id$
Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-09-15 16:03:37 UTC (rev 4553)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-09-15 17:38:57 UTC (rev 4554)
@@ -70,3 +70,38 @@
field->updateAdd(patch, *e_iter, _elemVector);
}
};
+
+void pylith::feassemble::Integrator::integrateLaplacianAction(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+{
+ const topology_type::patch_type patch = 0;
+ const Obj<topology_type>& topology = X->getTopology();
+ const Obj<topology_type::label_sequence>& elements = topology->heightStratum(patch, 0);
+ const topology_type::label_sequence::iterator end = elements->end();
+ value_type detJ;
+
+ for(topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != end; ++e_iter) {
+ computeElementGeometry(coordinates, *e_iter, _v0, _Jac, _invJac, detJ);
+ // Element integral
+ PetscMemzero(_elemVector, NUM_BASIS_FUNCTIONS*sizeof(value_type));
+ PetscMemzero(_elemMatrix, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(value_type));
+ for(int q = 0; q < NUM_QUADRATURE_POINTS; q++) {
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ _testWork[0] = _invJac[0]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*2+0] + _invJac[2]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*2+1];
+ _testWork[1] = _invJac[1]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*2+0] + _invJac[3]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*2+1];
+ for(int j = 0; j < NUM_BASIS_FUNCTIONS; j++) {
+ _basisWork[0] = _invJac[0]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*2+0] + _invJac[2]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*2+1];
+ _basisWork[1] = _invJac[1]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*2+0] + _invJac[3]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*2+1];
+ _elemMatrix[i*NUM_BASIS_FUNCTIONS+j] += (_testWork[0]*_basisWork[0] + _testWork[1]*_basisWork[1])*_weights[q]*detJ;
+ }
+ }
+ }
+ // Assembly
+ const value_type *xWork = X->restrict(patch, *e_iter);
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ for(int j = 0; j < NUM_BASIS_FUNCTIONS; j++) {
+ _elemVector[i] += _elemMatrix[i*NUM_BASIS_FUNCTIONS+j]*xWork[j];
+ }
+ }
+ F->updateAdd(patch, *e_iter, _elemVector);
+ }
+};
Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh 2006-09-15 16:03:37 UTC (rev 4553)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh 2006-09-15 17:38:57 UTC (rev 4554)
@@ -18,21 +18,29 @@
value_type *_points;
value_type *_weights;
value_type *_basis;
+ value_type *_basisDer;
value_type *_v0;
value_type *_Jac;
value_type *_invJac;
value_type *_realCoords;
value_type *_elemVector;
+ value_type *_elemMatrix;
+ value_type *_testWork;
+ value_type *_basisWork;
public:
- Integrator(int dim, int numQuadPoints, value_type points[], value_type weights[], int numBasisFuncs, value_type basis[]) : NUM_QUADRATURE_POINTS(numQuadPoints), NUM_BASIS_FUNCTIONS(numBasisFuncs), _dim(dim) {
+ Integrator(int dim, int numQuadPoints, value_type points[], value_type weights[], int numBasisFuncs, value_type basis[], value_type basisDer[]) : NUM_QUADRATURE_POINTS(numQuadPoints), NUM_BASIS_FUNCTIONS(numBasisFuncs), _dim(dim) {
_points = new value_type[numQuadPoints*dim];
_weights = new value_type[numQuadPoints];
_basis = new value_type[numQuadPoints*numBasisFuncs];
+ _basisDer = new value_type[numQuadPoints*numBasisFuncs*dim];
_v0 = new value_type[dim];
_Jac = new value_type[dim*dim];
_invJac = new value_type[dim*dim];
_realCoords = new value_type[dim];
_elemVector = new value_type[numBasisFuncs];
+ _elemMatrix = new value_type[numBasisFuncs*numBasisFuncs];
+ _testWork = new value_type[numBasisFuncs];
+ _basisWork = new value_type[numBasisFuncs];
for(int q = 0; q < numQuadPoints; ++q) {
for(int d = 0; d < dim; ++d) {
@@ -41,6 +49,9 @@
_weights[q] = weights[q];
for(int b = 0; b < numBasisFuncs; ++b) {
_basis[q*numBasisFuncs+b] = basis[q*numBasisFuncs+b];
+ for(int d = 0; d < dim; ++d) {
+ _basisDer[(q*numBasisFuncs+b)*dim+d] = basisDer[(q*numBasisFuncs+b)*dim+d];
+ }
}
}
};
@@ -53,10 +64,14 @@
delete [] _invJac;
delete [] _realCoords;
delete [] _elemVector;
+ delete [] _elemMatrix;
+ delete [] _testWork;
+ delete [] _basisWork;
};
public:
void computeElementGeometry(const Obj<section_type>&, const point_type&, value_type [], value_type [], value_type [], value_type&);
void integrateFunction(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(value_type []));
+ void integrateLaplacianAction(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
};
}
}
Modified: short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc 2006-09-15 16:03:37 UTC (rev 4553)
+++ short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc 2006-09-15 17:38:57 UTC (rev 4554)
@@ -10,14 +10,14 @@
// ======================================================================
//
-#include "MeshIOAscii.hh"
+#include "pylith/meshio/MeshIOAscii.hh"
#include "elemVector.hh"
#include <iostream> // USES std::cerr
#include "Integration.hh"
// ----------------------------------------------------------------------
-ALE::Mesh::section_type::value_type f(ALE::Mesh::section_type::value_type coords[])
+ALE::Mesh::section_type::value_type foo(ALE::Mesh::section_type::value_type coords[])
{
return 1.0;
}
@@ -38,20 +38,25 @@
try {
ALE::Obj<ALE::Mesh> mesh;
- pylith::meshIO::MeshIOAscii iohandler;
+ pylith::meshio::MeshIOAscii iohandler;
iohandler.filename(argv[1]);
iohandler.read(mesh, false);
const ALE::Mesh::topology_type::patch_type patch = 0;
- const Obj<ALE::Mesh::section_type>& field = mesh->getSection("field");
+ const Obj<ALE::Mesh::section_type>& X = mesh->getSection("X");
+ const Obj<ALE::Mesh::section_type>& F = mesh->getSection("F");
pylith::feassemble::Integrator integrator(mesh->getDimension(),
NUM_QUADRATURE_POINTS, points, weights,
- NUM_BASIS_FUNCTIONS, Basis);
+ NUM_BASIS_FUNCTIONS, Basis, BasisDerivatives);
- field->setFiberDimensionByDepth(patch, 0, 1);
- field->allocate();
- integrator.integrateFunction(field, mesh->getSection("coordinates"), f);
- field->view("Weak form of f");
+ X->setFiberDimensionByDepth(patch, 0, 1);
+ X->allocate();
+ F->setFiberDimensionByDepth(patch, 0, 1);
+ F->allocate();
+ integrator.integrateFunction(X, mesh->getSection("coordinates"), foo);
+ X->view("Weak form of foo");
+ integrator.integrateLaplacianAction(X, F, mesh->getSection("coordinates"));
+ F->view("Weak form of \Delta foo");
} catch(ALE::Exception e) {
int rank;
MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
More information about the cig-commits
mailing list