[cig-commits] r8376 - in short/3D/PyLith/trunk: . libsrc/feassemble
brad at geodynamics.org
brad at geodynamics.org
Tue Dec 4 11:32:26 PST 2007
Author: brad
Date: 2007-12-04 11:32:25 -0800 (Tue, 04 Dec 2007)
New Revision: 8376
Modified:
short/3D/PyLith/trunk/README
short/3D/PyLith/trunk/TODO
short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.hh
short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.cc
short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.hh
Log:
Updated README and TODO for 1.0.2 release.
Modified: short/3D/PyLith/trunk/README
===================================================================
--- short/3D/PyLith/trunk/README 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/README 2007-12-04 19:32:25 UTC (rev 8376)
@@ -38,11 +38,82 @@
and the generation of Green's functions to be used in inversions.
Many of these features should be available in Fall 2007.
+
======================================================================
+TIPS
+======================================================================
+
+ * For most crustal deformation problems, we recommend using the
+ Additive Schwartz preconditioner via the following PETSc
+ settings:
+
+ - Command line arguments
+
+ --petsc.pc_type=asm
+ --petsc.ksp_max_it=100
+ --petsc.ksp_gmres_restart=50
+ --petsc.ksp_rtol=1.0e-08
+
+ - pylithapp.cfg (or your other favorite .cfg file)
+
+ [pylithapp.petsc]
+ pc_type = asm
+ ksp_max_it = 200
+ ksp_gmres_restart = 50
+ ksp_rtol = 1.0e-08
+
+ * If the solve takes more than a few hundred iterations for a
+ large problem (use the --petsc.ksp_monitor=1 and
+ --petsc.ksp_view=1 to see diagnostic information for the
+ solver), this is an indication that something is wrong. Either
+ the preconditioner is inappropriate for the type of problem you
+ are solving or there is an error in the problem setup.
+
+
+======================================================================
RELEASE NOTES
======================================================================
----------------------------------------------------------------------
+Version 1.0.2
+----------------------------------------------------------------------
+
+ * Performance optimizations have significantly reduced runtime and
+ memory use relative to version 1.0.1. The default quadrature order
+ for tetrahedral cells is now 1, which is appropriate for the
+ default basis functions.
+
+ * Added checks to verify the compatibility of quadrature scheme for solid
+ and cohesive cells.
+
+ * Bug fixes
+
+ - In some cases, cohesive cells were not inserted into the
+ finite-element mesh properly. The cells mixed together vertices
+ from the different sides of the fault. A more efficient
+ procedure for creating cohesive cells fixed this problem.
+
+ - Cell adjacency graph was created incorrectly which resulted in
+ a poor quality of partitioning among processors.
+
+ - VTK output for meshes with N faults included cohesive cells
+ for N-1 faults. Since VTK output does not understand cohesive
+ cells, we now remove all cohesive cells from the VTK output.
+
+ - Using the SimpleDB in Spatialdata from Python limited
+ interpolation to the "linear" scheme instead of allowing use of
+ the "nearest" scheme. Setting the SimpleDB property to "nearest"
+ and "linear" now works as expected.
+
+ - The reader for Spatialdata coordinate systems information did not
+ correctly putback characters in the input stream, resulting in
+ reading errors. The putback routines were fixed.
+
+ - Fault "up" and "normal" directions remained as string arrays
+ when passed to the module, instead of being converted to float
+ arrays.
+
+----------------------------------------------------------------------
Version 1.0.1
----------------------------------------------------------------------
@@ -87,7 +158,7 @@
to fault slip, which increase the accuracy of displacement
fields near faults and facilitate visualization of fault slip.
- - Usage of cohesive elements will facilitiate the upcoming
+ - Usage of cohesive elements will facilitate the upcoming
addition of fault constitutive relations, where fault slip
occurs in response to prescribed physics.
@@ -105,6 +176,7 @@
- Generalized Maxwell and Power-law Maxwell viscoelastic models
+
======================================================================
KNOWN ISSUES
======================================================================
Modified: short/3D/PyLith/trunk/TODO
===================================================================
--- short/3D/PyLith/trunk/TODO 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/TODO 2007-12-04 19:32:25 UTC (rev 8376)
@@ -2,48 +2,18 @@
CURRENT ISSUES
======================================================================
-Release 1.0.2
-
- 1. Fix memory imbalance associated with distribution.
-
- 2. Fix problems with buildbot on darwin? (is binary ok?)
-
- 3. Add check to make sure every material in mesh has a material model.
-
- Add check for overlapping of material ids for bulk and cohesive
- cells.
-
Release 1.1
- 1. EventLogging
-
+ 1. Finish Neumann BC [CHARLES]
+
+ 2. EventLogging [BRAD]
Check constraints/integrators.
Fault::adjustTopology()
High level events
MeshGenerator
Problem
- 2. Optimization of restrict()/update()
-
- a. Fault
- b. Absorbing BC
- c. Neumann BC
-
- 3. Nonisoparametric cells
- Dirichlet and fault constraints?
- Fix Quadrature2DLinear.py (wrong quad pts and wts)?
- Cleanup CellGeometry
- Allow multiple locations for ptRefToGlobal() and jacobian() for
- better efficiency
- Switch from double_array to double* (used by faults)
-
- 4. Finish Neumann BC.
-
- 5. Finish velocity Dirichlet BC.
-
- 6. Generalized Maxwell viscoelastic model
-
- 7. Reimplement SolutionIO.
+ 3. Reimplement SolutionIO [BRAD & MATT]
a. Reimplement SolutionIOVTK
b. Implement SolutionIOHDF5
@@ -118,6 +88,15 @@
Solution field (time history)
+ 4. Optimization of restrict()/update() [BRAD]
+ a. Fault
+ b. Absorbing BC
+ c. Neumann BC
+
+ 5. Velocity Dirichlet BC [BRAD]
+
+ 6. Generalized Maxwell viscoelastic model [CHARLES]
+
General
1. Need to add explanation of output and output parameters to
@@ -125,8 +104,13 @@
2. Add dependency diagram to manual.
- 3. Experiment with different preconditioners and tabulate results
+ 3. Fix memory imbalance associated with distribution.
+ 4. Add check to make sure every material in mesh has a material model.
+
+ Add check for overlapping of material ids for bulk and cohesive
+ cells.
+
Release 1.2
1. Dynamic fault interface conditions.
@@ -135,8 +119,16 @@
3. Nonlinear solver
- 4. Replace use of Pyrexembed with SWIG
+ 4. Nonisoparametric cells
+ C++ unit tests for CellGeometry refPtsToGlobal() and jacobian()
+ Update quadrature to use CellGeometry refPtsToGlobal() and jacobian()
+ Remove CellGeometry jacobian(double_array)
+
+ Dirichlet and fault constraints?
+
+ 5. Replace use of Pyrexembed with SWIG
+
======================================================================
KNOWN DEFICIENCIES
======================================================================
Modified: short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/CellGeometry.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -101,16 +101,18 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
virtual
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const = 0;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const = 0;
/** Compute Jacobian at location in cell.
*
@@ -130,15 +132,17 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
virtual
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const = 0;
+ const double* ptsRef,
+ const int dim,
+ const int npts) const = 0;
/** Compute orientation of cell at location.
*
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -65,13 +65,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryHex3D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryHex3D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -108,13 +109,18 @@
const double y7 = vertices[22];
const double z7 = vertices[23];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- const double p1 = 0.5*(1.0+coordsRef[1]);
- const double p2 = 0.5*(1.0+coordsRef[2]);
- assert(0 <= p0 && p0 <= 1.0);
- assert(0 <= p1 && p1 <= 1.0);
- assert(0 <= p2 && p2 <= 1.0);
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
+ const double h_1 = z1 - z0;
+ const double f_3 = x3 - x0;
+ const double g_3 = y3 - y0;
+ const double h_3 = z3 - z0;
+
+ const double f_4 = x4 - x0;
+ const double g_4 = y4 - y0;
+ const double h_4 = z4 - z0;
+
const double f_01 = x2 - x1 - x3 + x0;
const double g_01 = y2 - y1 - y3 + y0;
const double h_01 = z2 - z1 - z3 + z0;
@@ -131,15 +137,23 @@
const double g_012 = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
const double h_012 = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
- coordsGlobal[0] = x0 + (x1-x0) * p0 + (x3-x0) * p1 + (x4-x0) * p2
- + f_01*p0*p1 + f_12*p1*p2 + f_02*p0*p2 + f_012*p0*p1*p2;
- coordsGlobal[1] = y0 + (y1-y0) * p0 + (y3-y0) * p1 + (y4-y0) * p2
- + g_01*p0*p1 + g_12*p1*p2 + g_02*p0*p2 + g_012*p0*p1*p2;
- coordsGlobal[2] = z0 + (z1-z0) * p0 + (z3-z0) * p1 + (z4-z0) * p2
- + h_01*p0*p1 + h_12*p1*p2 + h_02*p0*p2 + h_012*p0*p1*p2;
+ for (int i=0, iR=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p2 = 0.5 * (1.0 + ptsRef[iR++]);
+ assert(0 <= p0 && p0 <= 1.0);
+ assert(0 <= p1 && p1 <= 1.0);
+ assert(0 <= p2 && p2 <= 1.0);
+ ptsGlobal[iG++] = x0 + f_1*p0 + f_3*p1 + f_4*p2
+ + f_01*p0*p1 + f_12*p1*p2 + f_02*p0*p2 + f_012*p0*p1*p2;
+ ptsGlobal[iG++] = y0 + g_1*p0 + g_3*p1 + g_4*p2 + g_01*p0*p1
+ + g_01*p0*p1 + g_12*p1*p2 + g_02*p0*p2 + g_012*p0*p1*p2;
+ ptsGlobal[iG++] = z0 + h_1*p0 + h_3*p1 + h_4*p2 + h_01*p0*p1
+ + h_01*p0*p1 + h_12*p1*p2 + h_02*p0*p2 + h_012*p0*p1*p2;
+ } // for
- PetscLogFlopsNoCheck(120);
-} // coordsRefToGlobal
+ PetscLogFlopsNoCheck(57 + npts*57);
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
@@ -240,13 +254,14 @@
pylith::feassemble::GeometryHex3D::jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const
+ const double* ptsRef,
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
assert(0 != vertices);
- assert(0 != location);
+ assert(0 != ptsRef);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -282,47 +297,66 @@
const double y7 = vertices[22];
const double z7 = vertices[23];
- const double x = 0.5 * (location[0] + 1.0);
- const double y = 0.5 * (location[1] + 1.0);
- const double z = 0.5 * (location[2] + 1.0);
- assert(-1.0 <= x && x <= 1.0);
- assert(-1.0 <= y && y <= 1.0);
- assert(-1.0 <= z && z <= 1.0);
+ const double f_1 = (x1 - x0) / 2.0;
+ const double g_1 = (y1 - y0) / 2.0;
+ const double h_1 = (z1 - z0) / 2.0;
- const double f_xy = x2 - x1 - x3 + x0;
- const double g_xy = y2 - y1 - y3 + y0;
- const double h_xy = z2 - z1 - z3 + z0;
+ const double f_3 = (x3 - x0) / 2.0;
+ const double g_3 = (y3 - y0) / 2.0;
+ const double h_3 = (z3 - z0) / 2.0;
+
+ const double f_4 = (x4 - x0) / 2.0;
+ const double g_4 = (y4 - y0) / 2.0;
+ const double h_4 = (z4 - z0) / 2.0;
- const double f_yz = x7 - x3 - x4 + x0;
- const double g_yz = y7 - y3 - y4 + y0;
- const double h_yz = z7 - z3 - z4 + z0;
+ const double f_01 = (x2 - x1 - x3 + x0) / 2.0;
+ const double g_01 = (y2 - y1 - y3 + y0) / 2.0;
+ const double h_01 = (z2 - z1 - z3 + z0) / 2.0;
- const double f_xz = x5 - x1 - x4 + x0;
- const double g_xz = y5 - y1 - y4 + y0;
- const double h_xz = z5 - z1 - z4 + z0;
+ const double f_12 = (x7 - x3 - x4 + x0) / 2.0;
+ const double g_12 = (y7 - y3 - y4 + y0) / 2.0;
+ const double h_12 = (z7 - z3 - z4 + z0) / 2.0;
- const double f_xyz = x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7;
- const double g_xyz = y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7;
- const double h_xyz = z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7;
+ const double f_02 = (x5 - x1 - x4 + x0) / 2.0;
+ const double g_02 = (y5 - y1 - y4 + y0) / 2.0;
+ const double h_02 = (z5 - z1 - z4 + z0) / 2.0;
- jacobian[0] = (x1 - x0 + f_xy*y + f_xz*z + f_xyz*y*z) / 2.0;
- jacobian[1] = (x3 - x0 + f_xy*x + f_yz*z + f_xyz*x*z) / 2.0;
- jacobian[2] = (x4 - x0 + f_yz*y + f_xz*x + f_xyz*x*y) / 2.0;
+ const double f_012 = (x6 - x0 + x1 - x2 + x3 + x4 - x5 - x7) / 2.0;
+ const double g_012 = (y6 - y0 + y1 - y2 + y3 + y4 - y5 - y7) / 2.0;
+ const double h_012 = (z6 - z0 + z1 - z2 + z3 + z4 - z5 - z7) / 2.0;
- jacobian[3] = (y1 - y0 + g_xy*y + g_xz*z + g_xyz*y*z) / 2.0;
- jacobian[4] = (y3 - y0 + g_xy*x + g_yz*z + g_xyz*x*z) / 2.0;
- jacobian[5] = (y4 - y0 + g_yz*y + g_xz*x + g_xyz*x*y) / 2.0;
+ for (int i=0, iR=0, iJ=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p2 = 0.5 * (1.0 + ptsRef[iR++]);
+ assert(0 <= p0 && p0 <= 1.0);
+ assert(0 <= p1 && p1 <= 1.0);
+ assert(0 <= p2 && p2 <= 1.0);
+ const double j0 = f_1 + f_01*p1 + f_02*p2 + f_012*p1*p2;
+ const double j1 = f_3 + f_01*p0 + f_12*p2 + f_012*p0*p2;
+ const double j2 = f_4 + f_12*p1 + f_02*p0 + f_012*p0*p1;
+
+ const double j3 = g_1 + g_01*p1 + g_02*p2 + g_012*p1*p2;
+ const double j4 = g_3 + g_01*p0 + g_12*p2 + g_012*p0*p2;
+ const double j5 = g_4 + g_12*p1 + g_02*p0 + g_012*p0*p1;
+
+ const double j6 = h_1 + h_01*p1 + h_02*p2 + h_012*p1*p2;
+ const double j7 = h_3 + h_01*p0 + h_12*p2 + h_012*p0*p2;
+ const double j8 = h_4 + h_12*p1 + h_02*p0 + h_012*p0*p1;
- jacobian[6] = (z1 - z0 + h_xy*y + h_xz*z + h_xyz*y*z) / 2.0;
- jacobian[7] = (z3 - z0 + h_xy*x + h_yz*z + h_xyz*x*z) / 2.0;
- jacobian[8] = (z4 - z0 + h_yz*y + h_xz*x + h_xyz*x*y) / 2.0;
+ jacobian[iJ++] = j0;
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ jacobian[iJ++] = j4;
+ jacobian[iJ++] = j5;
+ jacobian[iJ++] = j6;
+ jacobian[iJ++] = j7;
+ jacobian[iJ++] = j8;
+ det[i] = j0*(j4*j8 - j5*j7) - j1*(j3*j8 - j5*j6) + j2*(j3*j7 - j4*j6);
+ } // for
- *det =
- jacobian[0]*(jacobian[4]*jacobian[8] - jacobian[5]*jacobian[7]) -
- jacobian[1]*(jacobian[3]*jacobian[8] - jacobian[5]*jacobian[6]) +
- jacobian[2]*(jacobian[3]*jacobian[7] - jacobian[4]*jacobian[6]);
-
- PetscLogFlopsNoCheck(152);
+ PetscLogFlopsNoCheck(78 + npts*69);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryHex3D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -67,15 +67,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -94,14 +96,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -59,13 +59,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryLine1D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryLine1D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(1 == dim);
assert(spaceDim() == dim);
@@ -73,11 +74,13 @@
const double x0 = vertices[0];
const double x1 = vertices[1];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- coordsGlobal[0] = x0 + (x1-x0) * p0;
+ for (int i=0; i < npts; ++i) {
+ const double p0 = 0.5*(1.0+ptsRef[i]);
+ ptsGlobal[i] = x0 + (x1-x0) * p0;
+ } // for
- PetscLogFlopsNoCheck(5);
-} // coordsRefToGlobal
+ PetscLogFlopsNoCheck(5*npts);
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
@@ -107,7 +110,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -119,8 +123,12 @@
const double x0 = vertices[0];
const double x1 = vertices[1];
- jacobian[0] = (x1 - x0)/2.0;
- *det = jacobian[0];
+ const double j = (x1 - x0)/2.0;
+ for (int i=0; i < npts; ++i) {
+ jacobian[i] = j;
+ det[i] = j;
+ } // for
+
PetscLogFlopsNoCheck(2);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine1D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -55,15 +55,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -82,14 +84,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -59,13 +59,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryLine2D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryLine2D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -76,20 +77,25 @@
const double x1 = vertices[2];
const double y1 = vertices[3];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- coordsGlobal[0] = x0 + (x1-x0) * p0;
- coordsGlobal[1] = y0 + (y1-y0) * p0;
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
- PetscLogFlopsNoCheck(8);
-} // coordsRefToGlobal
+ for (int i=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[i]);
+ ptsGlobal[iG++] = x0 + f_1 * p0;
+ ptsGlobal[iG++] = y0 + g_1 * p0;
+ } // for
+ PetscLogFlopsNoCheck(2 + npts*6);
+} // ptsRefToGlobal
+
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryLine2D::jacobian(double_array* jacobian,
- double* det,
- const double_array& vertices,
- const double_array& location) const
+ double* det,
+ const double_array& vertices,
+ const double_array& location) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -118,7 +124,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -133,10 +140,16 @@
const double x1 = vertices[2];
const double y1 = vertices[3];
- jacobian[0] = (x1 - x0)/2.0;
- jacobian[1] = (y1 - y0)/2.0;
- *det = sqrt(pow(jacobian[0], 2) +
- pow(jacobian[1], 2));
+ const double j1 = (x1 - x0) / 2.0;
+ const double j2 = (y1 - y0) / 2.0;
+ const double jdet = sqrt(j1*j1 + j2*j2);
+
+ for (int i=0, iJ=0; i < npts; ++i) {
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ det[i] = jdet;
+ } // for
+
PetscLogFlopsNoCheck(8);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine2D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -55,15 +55,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -82,14 +84,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -59,13 +59,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryLine3D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryLine3D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -78,14 +79,20 @@
const double y1 = vertices[4];
const double z1 = vertices[5];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- coordsGlobal[0] = x0 + (x1-x0) * p0;
- coordsGlobal[1] = y0 + (y1-y0) * p0;
- coordsGlobal[2] = z0 + (z1-z0) * p0;
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
+ const double h_1 = z1 - z0;
- PetscLogFlopsNoCheck(11);
-} // coordsRefToGlobal
+ for (int i=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[i]);
+ ptsGlobal[iG++] = x0 + f_1 * p0;
+ ptsGlobal[iG++] = y0 + g_1 * p0;
+ ptsGlobal[iG++] = z0 + h_1 * p0;
+ } // for
+ PetscLogFlopsNoCheck(3 + npts*8);
+} // ptsRefToGlobal
+
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
@@ -124,7 +131,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -141,12 +149,18 @@
const double y1 = vertices[4];
const double z1 = vertices[5];
- jacobian[0] = (x1 - x0)/2.0;
- jacobian[1] = (y1 - y0)/2.0;
- jacobian[2] = (z1 - z0)/2.0;
- *det = sqrt(pow(jacobian[0], 2) +
- pow(jacobian[1], 2) +
- pow(jacobian[2], 2));
+ const double j1 = (x1 - x0) / 2.0;
+ const double j2 = (y1 - y0) / 2.0;
+ const double j3 = (z1 - z0) / 2.0;
+ const double jdet = sqrt(j1*j1 + j2*j2 + j3*j3);
+
+ for (int i=0, iJ=0; i < npts; ++i) {
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ det[i] = jdet;
+ } // for
+
PetscLogFlopsNoCheck(12);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryLine3D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -55,15 +55,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -82,14 +84,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -52,27 +52,29 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryPoint1D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryPoint1D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(1 == dim);
assert(spaceDim() == dim);
- coordsGlobal[0] = vertices[0];
-} // coordsRefToGlobal
+ for (int i=0; i < npts; ++i)
+ ptsGlobal[i] = vertices[i];
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryPoint1D::jacobian(double_array* jacobian,
- double* det,
- const double_array& vertices,
- const double_array& location) const
+ double* det,
+ const double_array& vertices,
+ const double_array& location) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -89,18 +91,21 @@
pylith::feassemble::GeometryPoint1D::jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const
+ const double* ptsRef,
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
assert(0 != vertices);
- assert(0 != location);
+ assert(0 != ptsRef);
assert(1 == dim);
assert(spaceDim() == dim);
- jacobian[0] = 1.0;
- *det = 1.0;
+ for (int i=0; i < npts; ++i) {
+ jacobian[i] = 1.0;
+ det[i] = 1.0;
+ } // for
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint1D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -55,15 +55,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -82,14 +84,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -52,20 +52,22 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryPoint2D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryPoint2D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(2 == dim);
assert(spaceDim() == dim);
- coordsGlobal[0] = vertices[0];
- coordsGlobal[1] = vertices[1];
-} // coordsRefToGlobal
+ const int size = npts*dim;
+ for (int i=0; i < size; ++i)
+ ptsGlobal[i] = vertices[i];
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
@@ -91,7 +93,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -100,8 +103,10 @@
assert(2 == dim);
assert(spaceDim() == dim);
- jacobian[0] = 1.0;
- *det = 1.0;
+ for (int i=0; i < npts; ++i) {
+ jacobian[i] = 1.0;
+ det[i] = 1.0;
+ } // for
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint2D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -55,15 +55,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -82,14 +84,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -52,21 +52,22 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryPoint3D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryPoint3D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(3 == dim);
assert(spaceDim() == dim);
- coordsGlobal[0] = vertices[0];
- coordsGlobal[1] = vertices[1];
- coordsGlobal[2] = vertices[2];
-} // coordsRefToGlobal
+ const int size = npts*dim;
+ for (int i=0; i < size; ++i)
+ ptsGlobal[i] = vertices[i];
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
@@ -92,7 +93,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -101,8 +103,10 @@
assert(3 == dim);
assert(spaceDim() == dim);
- jacobian[0] = 1.0;
- *det = 1.0;
+ for (int i=0; i < npts; ++i) {
+ jacobian[i] = 1.0;
+ det[i] = 1.0;
+ } // for
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryPoint3D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -55,15 +55,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -82,14 +84,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -61,13 +61,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryQuad2D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryQuad2D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -84,16 +85,25 @@
const double x3 = vertices[6];
const double y3 = vertices[7];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- const double p1 = 0.5*(1.0+coordsRef[1]);
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
+
+ const double f_3 = x3 - x0;
+ const double g_3 = y3 - y0;
+
const double f_01 = x2 - x1 - x3 + x0;
const double g_01 = y2 - y1 - y3 + y0;
- coordsGlobal[0] = x0 + (x1-x0) * p0 + (x2-x0) * p1 + f_01 * p0 * p1;
- coordsGlobal[1] = y0 + (y1-y0) * p0 + (y2-y0) * p1 + g_01 * p0 * p1;
- PetscLogFlopsNoCheck(28);
-} // coordsRefToGlobal
+ for (int i=0, iR=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ ptsGlobal[iG++] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
+ ptsGlobal[iG++] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
+ } // for
+ PetscLogFlopsNoCheck(10 + npts*18);
+} // ptsRefToGlobal
+
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
@@ -147,13 +157,14 @@
pylith::feassemble::GeometryQuad2D::jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const
+ const double* ptsRef,
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
assert(0 != vertices);
- assert(0 != location);
+ assert(0 != ptsRef);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -169,23 +180,33 @@
const double x3 = vertices[6];
const double y3 = vertices[7];
- const double x = 0.5 * (location[0] + 1.0);
- const double y = 0.5 * (location[1] + 1.0);
- assert(0 <= x && x <= 1.0);
- assert(0 <= y && y <= 1.0);
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
- const double f_xy = x2 - x1 - x3 + x0;
- const double g_xy = y2 - y1 - y3 + y0;
+ const double f_3 = x3 - x0;
+ const double g_3 = y3 - y0;
- jacobian[0] = (x1 - x0 + f_xy*y) / 2.0;
- jacobian[1] = (x3 - x0 + f_xy*x) / 2.0;
- jacobian[2] = (y1 - y0 + g_xy*y) / 2.0;
- jacobian[3] = (y3 - y0 + g_xy*x) / 2.0;
+ const double f_01 = x2 - x1 - x3 + x0;
+ const double g_01 = y2 - y1 - y3 + y0;
- *det =
- jacobian[0]*jacobian[3] -
- jacobian[1]*jacobian[2];
- PetscLogFlopsNoCheck(29);
+ for (int i=0, iR=0, iJ=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ assert(0 <= p0 && p0 <= 1.0);
+ assert(0 <= p1 && p1 <= 1.0);
+ const double j0 = (f_1 + f_01 * p1) / 2.0;
+ const double j1 = (f_3 + f_01 * p0) / 2.0;
+ const double j2 = (g_1 + g_01 * p1) / 2.0;
+ const double j3 = (g_3 + g_01 * p0) / 2.0;
+
+ jacobian[iJ++] = j0;
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ det[i] = j0*j3 - j1*j2;
+ } // for
+
+ PetscLogFlopsNoCheck(10 + npts*19);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad2D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -68,15 +68,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -95,14 +97,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -61,13 +61,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryQuad3D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryQuad3D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -88,20 +89,30 @@
const double y3 = vertices[10];
const double z3 = vertices[11];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- const double p1 = 0.5*(1.0+coordsRef[1]);
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
+ const double h_1 = z1 - z0;
+
+ const double f_3 = x3 - x0;
+ const double g_3 = y3 - y0;
+ const double h_3 = z3 - z0;
+
const double f_01 = x2 - x1 - x3 + x0;
const double g_01 = y2 - y1 - y3 + y0;
const double h_01 = z2 - z1 - z3 + z0;
- assert(0 <= p0 && p0 <= 1.0);
- assert(0 <= p1 && p1 <= 1.0);
- coordsGlobal[0] = x0 + (x1-x0) * p0 + (x2-x0) * p1 + f_01 * p0 * p1;
- coordsGlobal[1] = y0 + (y1-y0) * p0 + (y2-y0) * p1 + g_01 * p0 * p1;
- coordsGlobal[2] = z0 + (z1-z0) * p0 + (z2-z0) * p1 + h_01 * p0 * p1;
+ for (int i=0, iR=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ assert(0 <= p0 && p0 <= 1.0);
+ assert(0 <= p1 && p1 <= 1.0);
+ ptsGlobal[iG++] = x0 + f_1 * p0 + f_3 * p1 + f_01 * p0 * p1;
+ ptsGlobal[iG++] = y0 + g_1 * p0 + g_3 * p1 + g_01 * p0 * p1;
+ ptsGlobal[iG++] = z0 + h_1 * p0 + h_3 * p1 + h_01 * p0 * p1;
+ } // for
- PetscLogFlopsNoCheck(40);
-} // coordsRefToGlobal
+ PetscLogFlopsNoCheck(15 + npts*25);
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
@@ -176,13 +187,14 @@
pylith::feassemble::GeometryQuad3D::jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const
+ const double* ptsRef,
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
assert(0 != vertices);
- assert(0 != location);
+ assert(0 != ptsRef);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -202,40 +214,44 @@
const double y3 = vertices[10];
const double z3 = vertices[11];
- const double x = 0.5 * (location[0] + 1.0);
- const double y = 0.5 * (location[1] + 1.0);
- assert(0 <= x && x <= 1.0);
- assert(0 <= y && y <= 1.0);
+ const double f_1 = (x1 - x0) / 2.0;
+ const double g_1 = (y1 - y0) / 2.0;
+ const double h_1 = (z1 - z0) / 2.0;
- const double f_xy = x2 - x1 - x3 + x0;
- const double g_xy = y2 - y1 - y3 + y0;
- const double h_xy = z2 - z1 - z3 + z0;
+ const double f_3 = (x3 - x0) / 2.0;
+ const double g_3 = (y3 - y0) / 2.0;
+ const double h_3 = (z3 - z0) / 2.0;
- jacobian[0] = (x1 - x0 + f_xy*y) / 2.0;
- jacobian[1] = (x3 - x0 + f_xy*x) / 2.0;
+ const double f_01 = (x2 - x1 - x3 + x0) / 2.0;
+ const double g_01 = (y2 - y1 - y3 + y0) / 2.0;
+ const double h_01 = (z2 - z1 - z3 + z0) / 2.0;
- jacobian[2] = (y1 - y0 + g_xy*y) / 2.0;
- jacobian[3] = (y3 - y0 + g_xy*x) / 2.0;
+ for (int i=0, iR=0, iJ=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ assert(0 <= p0 && p0 <= 1.0);
+ assert(0 <= p1 && p1 <= 1.0);
+ const double j0 = f_1 + f_01 * p1;
+ const double j1 = f_3 + f_01 * p0;
+ const double j2 = g_1 + g_01 * p1;
+ const double j3 = g_3 + g_01 * p0;
+ const double j4 = h_1 + h_01 * p1;
+ const double j5 = h_3 + h_01 * p0;
+ jacobian[iJ++] = j0;
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ jacobian[iJ++] = j4;
+ jacobian[iJ++] = j5;
- jacobian[4] = (z1 - z0 + h_xy*y) / 2.0;
- jacobian[5] = (z3 - z0 + h_xy*x) / 2.0;
+ const double jj00 = j0*j0 + j2*j2 + j4*j4;
+ const double jj10 = j0*j1 + j2*j3 + j4*j5;
+ const double jj01 = jj10;
+ const double jj11 = j1*j1 + j3*j3 + j5*j5;
+ det[i] = sqrt(jj00*jj11 - jj01*jj10);
+ } // for
- const double jj00 =
- jacobian[0]*jacobian[0] +
- jacobian[2]*jacobian[2] +
- jacobian[4]*jacobian[4];
- const double jj10 =
- jacobian[0]*jacobian[1] +
- jacobian[2]*jacobian[3] +
- jacobian[4]*jacobian[5];
- const double jj01 = jj10;
- const double jj11 =
- jacobian[1]*jacobian[1] +
- jacobian[3]*jacobian[3] +
- jacobian[5]*jacobian[5];
- *det = sqrt(jj00*jj11 - jj01*jj10);
-
- PetscLogFlopsNoCheck(50);
+ PetscLogFlopsNoCheck(28 + npts*32);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryQuad3D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -68,15 +68,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -95,14 +97,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -61,13 +61,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryTet3D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryTet3D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -88,27 +89,41 @@
const double y3 = vertices[10];
const double z3 = vertices[11];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- const double p1 = 0.5*(1.0+coordsRef[1]);
- const double p2 = 0.5*(1.0+coordsRef[2]);
- assert(0 <= p0 && p0 <= 1.0);
- assert(0 <= p1 && p1 <= 1.0);
- assert(0 <= p2 && p2 <= 1.0);
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
+ const double h_1 = z1 - z0;
- coordsGlobal[0] = x0 + (x1-x0) * p0 + (x2-x0) * p1 + (x3-x0) * p2;
- coordsGlobal[1] = y0 + (y1-y0) * p0 + (y2-y0) * p1 + (y3-y0) * p2;
- coordsGlobal[2] = z0 + (z1-z0) * p0 + (z2-z0) * p1 + (z3-z0) * p2;
+ const double f_2 = x2 - x0;
+ const double g_2 = y2 - y0;
+ const double h_2 = z2 - z0;
- PetscLogFlopsNoCheck(33);
-} // coordsRefToGlobal
+ const double f_3 = x3 - x0;
+ const double g_3 = y3 - y0;
+ const double h_3 = z3 - z0;
+ for (int i=0, iR=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p2 = 0.5 * (1.0 + ptsRef[iR++]);
+ assert(0 <= p0 && p0 <= 1.0);
+ assert(0 <= p1 && p1 <= 1.0);
+ assert(0 <= p2 && p2 <= 1.0);
+
+ ptsGlobal[iG++] = x0 + f_1 * p0 + f_2 * p1 + f_3 * p2;
+ ptsGlobal[iG++] = y0 + g_1 * p0 + g_2 * p1 + g_3 * p2;;
+ ptsGlobal[iG++] = z0 + h_1 * p0 + h_2 * p1 + h_3 * p2;
+ } // for
+
+ PetscLogFlopsNoCheck(9 + npts*24);
+} // ptsRefToGlobal
+
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryTet3D::jacobian(double_array* jacobian,
- double* det,
- const double_array& vertices,
- const double_array& location) const
+ double* det,
+ const double_array& vertices,
+ const double_array& location) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -159,13 +174,14 @@
pylith::feassemble::GeometryTet3D::jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const
+ const double* ptsRef,
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
assert(0 != vertices);
- assert(0 != location);
+ assert(0 != ptsRef);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -185,21 +201,33 @@
const double y3 = vertices[10];
const double z3 = vertices[11];
- jacobian[0] = (x1 - x0) / 2.0;
- jacobian[1] = (x2 - x0) / 2.0;
- jacobian[2] = (x3 - x0) / 2.0;
- jacobian[3] = (y1 - y0) / 2.0;
- jacobian[4] = (y2 - y0) / 2.0;
- jacobian[5] = (y3 - y0) / 2.0;
- jacobian[6] = (z1 - z0) / 2.0;
- jacobian[7] = (z2 - z0) / 2.0;
- jacobian[8] = (z3 - z0) / 2.0;
+ const double j0 = (x1 - x0) / 2.0;
+ const double j1 = (x2 - x0) / 2.0;
+ const double j2 = (x3 - x0) / 2.0;
- *det =
- jacobian[0]*(jacobian[4]*jacobian[8] - jacobian[5]*jacobian[7]) -
- jacobian[1]*(jacobian[3]*jacobian[8] - jacobian[5]*jacobian[6]) +
- jacobian[2]*(jacobian[3]*jacobian[7] - jacobian[4]*jacobian[6]);
+ const double j3 = (y1 - y0) / 2.0;
+ const double j4 = (y2 - y0) / 2.0;
+ const double j5 = (y3 - y0) / 2.0;
+ const double j6 = (z1 - z0) / 2.0;
+ const double j7 = (z2 - z0) / 2.0;
+ const double j8 = (z3 - z0) / 2.0;
+
+ const double jdet =
+ j0*(j4*j8 - j5*j7) -
+ j1*(j3*j8 - j5*j6) +
+ j2*(j3*j7 - j4*j6);
+
+ for (int i=0, iJ=0; i < npts; ++i) {
+ jacobian[iJ++] = j0;
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ jacobian[iJ++] = j4;
+ jacobian[iJ++] = j5;
+ det[i] = jdet;
+ } // for
+
PetscLogFlopsNoCheck(32);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryTet3D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -63,15 +63,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -90,14 +92,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -60,13 +60,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryTri2D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryTri2D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(2 == dim);
assert(spaceDim() == dim);
@@ -80,14 +81,22 @@
const double x2 = vertices[4];
const double y2 = vertices[5];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- const double p1 = 0.5*(1.0+coordsRef[1]);
- coordsGlobal[0] = x0 + (x1-x0) * p0 + (x2-x0) * p1;
- coordsGlobal[1] = y0 + (y1-y0) * p0 + (y2-y0) * p1;
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
- PetscLogFlopsNoCheck(16);
-} // coordsRefToGlobal
+ const double f_2 = x2 - x0;
+ const double g_2 = y2 - y0;
+ for (int i=0, iR=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ ptsGlobal[iG++] = x0 + f_1 * p0 + f_2 * p1;
+ ptsGlobal[iG++] = y0 + g_1 * p0 + g_2 * p1;
+ } // for
+
+ PetscLogFlopsNoCheck(4 + npts*12);
+} // ptsRefToGlobal
+
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
@@ -129,7 +138,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -147,15 +157,23 @@
const double x2 = vertices[4];
const double y2 = vertices[5];
- jacobian[0] = (x1 - x0) / 2.0;
- jacobian[1] = (x2 - x0) / 2.0;
- jacobian[2] = (y1 - y0) / 2.0;
- jacobian[3] = (y2 - y0) / 2.0;
- *det =
+ const double j1 = (x1 - x0) / 2.0;
+ const double j2 = (x2 - x0) / 2.0;
+ const double j3 = (y1 - y0) / 2.0;
+ const double j4 = (y2 - y0) / 2.0;
+ const double jdet =
jacobian[0]*jacobian[3] -
jacobian[1]*jacobian[2];
+ for (int i=0, iJ=0; i < npts; ++i) {
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ jacobian[iJ++] = j4;
+ det[i] = jdet;
+ } // for
+
PetscLogFlopsNoCheck(11);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri2D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -62,15 +62,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -89,14 +91,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.cc 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.cc 2007-12-04 19:32:25 UTC (rev 8376)
@@ -60,13 +60,14 @@
// ----------------------------------------------------------------------
// Transform coordinates in reference cell to global coordinates.
void
-pylith::feassemble::GeometryTri3D::coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const
-{ // coordsRefToGlobal
- assert(0 != coordsGlobal);
- assert(0 != coordsRef);
+pylith::feassemble::GeometryTri3D::ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts) const
+{ // ptsRefToGlobal
+ assert(0 != ptsGlobal);
+ assert(0 != ptsRef);
assert(0 != vertices);
assert(3 == dim);
assert(spaceDim() == dim);
@@ -83,22 +84,32 @@
const double y2 = vertices[7];
const double z2 = vertices[8];
- const double p0 = 0.5*(1.0+coordsRef[0]);
- const double p1 = 0.5*(1.0+coordsRef[1]);
- coordsGlobal[0] = x0 + (x1-x0) * p0 + (x2-x0) * p1;
- coordsGlobal[1] = y0 + (y1-y0) * p0 + (y2-y0) * p1;
- coordsGlobal[2] = z0 + (z1-z0) * p0 + (z2-z0) * p1;
+ const double f_1 = x1 - x0;
+ const double g_1 = y1 - y0;
+ const double h_1 = z1 - z0;
+ const double f_2 = x2 - x0;
+ const double g_2 = y2 - y0;
+ const double h_2 = z2 - z0;
+
+ for (int i=0, iR=0, iG=0; i < npts; ++i) {
+ const double p0 = 0.5 * (1.0 + ptsRef[iR++]);
+ const double p1 = 0.5 * (1.0 + ptsRef[iR++]);
+ ptsGlobal[iG++] = x0 + f_1 * p0 + f_2 * p1;
+ ptsGlobal[iG++] = y0 + g_1 * p0 + g_2 * p1;
+ ptsGlobal[iG++] = z0 + h_1 * p0 + h_2 * p1;
+ } // for
+
PetscLogFlopsNoCheck(22);
-} // coordsRefToGlobal
+} // ptsRefToGlobal
// ----------------------------------------------------------------------
// Compute Jacobian at location in cell.
void
pylith::feassemble::GeometryTri3D::jacobian(double_array* jacobian,
- double* det,
- const double_array& vertices,
- const double_array& location) const
+ double* det,
+ const double_array& vertices,
+ const double_array& location) const
{ // jacobian
assert(0 != jacobian);
@@ -150,7 +161,8 @@
double* det,
const double* vertices,
const double* location,
- const int dim) const
+ const int dim,
+ const int npts) const
{ // jacobian
assert(0 != jacobian);
assert(0 != det);
@@ -171,30 +183,31 @@
const double y2 = vertices[7];
const double z2 = vertices[8];
- jacobian[0] = (x1 - x0) / 2.0;
- jacobian[1] = (x2 - x0) / 2.0;
+ const double j0 = (x1 - x0) / 2.0;
+ const double j1 = (x2 - x0) / 2.0;
- jacobian[2] = (y1 - y0) / 2.0;
- jacobian[3] = (y2 - y0) / 2.0;
+ const double j2 = (y1 - y0) / 2.0;
+ const double j3 = (y2 - y0) / 2.0;
- jacobian[4] = (z1 - z0) / 2.0;
- jacobian[5] = (z2 - z0) / 2.0;
+ const double j4 = (z1 - z0) / 2.0;
+ const double j5 = (z2 - z0) / 2.0;
- const double jj00 =
- jacobian[0]*jacobian[0] +
- jacobian[2]*jacobian[2] +
- jacobian[4]*jacobian[4];
- const double jj10 =
- jacobian[0]*jacobian[1] +
- jacobian[2]*jacobian[3] +
- jacobian[4]*jacobian[5];
+ const double jj00 = j0*j0 + j2*j2 + j4*j4;
+ const double jj10 = j0*j1 + j2*j3 + j4*j5;
const double jj01 = jj10;
- const double jj11 =
- jacobian[1]*jacobian[1] +
- jacobian[3]*jacobian[3] +
- jacobian[5]*jacobian[5];
- *det = sqrt(jj00*jj11 - jj01*jj10);
+ const double jj11 = j1*j1 + j3*j3 + j5*j5;
+ const double jdet = sqrt(jj00*jj11 - jj01*jj10);
+ for (int i=0, iJ=0; i < npts; ++i) {
+ jacobian[iJ++] = j0;
+ jacobian[iJ++] = j1;
+ jacobian[iJ++] = j2;
+ jacobian[iJ++] = j3;
+ jacobian[iJ++] = j4;
+ jacobian[iJ++] = j5;
+ det[i] = jdet;
+ } // for
+
PetscLogFlopsNoCheck(31);
} // jacobian
Modified: short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.hh 2007-12-04 06:53:11 UTC (rev 8375)
+++ short/3D/PyLith/trunk/libsrc/feassemble/GeometryTri3D.hh 2007-12-04 19:32:25 UTC (rev 8376)
@@ -62,15 +62,17 @@
/** Transform coordinates in reference cell to global coordinates.
*
- * @param coordsGlobal Coordinates in global coordinate system.
- * @param coordsRef Coordinates in reference cell.
+ * @param ptsGlobal Array of points in global coordinate system.
+ * @param ptsRef Array of points in reference cell.
* @param vertices Array of cell vertices in global coordinates.
* @param dim Dimension of global coordinate system.
+ * @param npts Number of points to transform.
*/
- void coordsRefToGlobal(double* coordsGlobal,
- const double* coordsRef,
- const double* vertices,
- const int dim) const;
+ void ptsRefToGlobal(double* ptsGlobal,
+ const double* ptsRef,
+ const double* vertices,
+ const int dim,
+ const int npts =1) const;
/** Compute Jacobian at location in cell.
*
@@ -89,14 +91,16 @@
* @param jacobian Jacobian at location.
* @param det Determinant of Jacobian at location.
* @param vertices Coordinates of vertices of cell.
- * @param location Location in reference cell at which to compute Jacobian.
+ * @param ptsRef Points in reference cell at which to compute Jacobian.
* @param dim Dimension of coordinate system.
+ * @param npts Number of points to transform.
*/
void jacobian(double* jacobian,
double* det,
const double* vertices,
- const double* location,
- const int dim) const;
+ const double* ptsRef,
+ const int dim,
+ const int npts =1) const;
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
More information about the cig-commits
mailing list