[cig-commits] r15073 - in short/3D/PyLith/trunk: examples/3d/hex8 libsrc/meshio libsrc/topology modulesrc/topology pylith/meshio pylith/perf pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Wed May 27 19:03:26 PDT 2009
Author: brad
Date: 2009-05-27 19:03:25 -0700 (Wed, 27 May 2009)
New Revision: 15073
Modified:
short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc
short/3D/PyLith/trunk/libsrc/topology/Field.cc
short/3D/PyLith/trunk/libsrc/topology/Field.hh
short/3D/PyLith/trunk/modulesrc/topology/Field.i
short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubMesh.py
short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubSubMesh.py
short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py
short/3D/PyLith/trunk/pylith/problems/Explicit.py
short/3D/PyLith/trunk/pylith/problems/Formulation.py
short/3D/PyLith/trunk/pylith/problems/Implicit.py
Log:
Fixed bug related to dimensionalization of displacement field. Create field using chart. Fixed some more nondimensionalization bugs.
Modified: short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg
===================================================================
--- short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/examples/3d/hex8/pylithapp.cfg 2009-05-28 02:03:25 UTC (rev 15073)
@@ -37,9 +37,6 @@
# ----------------------------------------------------------------------
# problem
# ----------------------------------------------------------------------
-[pylithapp.timedependent]
-#normalizer = spatialdata.units.NondimElasticQuasistatic
-
[pylithapp.timedependent.formulation.time_step]
# Define the total time for the simulation and the default time step size.
total_time = 0.0*s ; total time of simulation
Modified: short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/libsrc/meshio/CellFilterAvg.cc 2009-05-28 02:03:25 UTC (rev 15073)
@@ -95,7 +95,7 @@
assert(0 != _fieldAvg);
_fieldAvg->newSection(sectionIn->getChart(), fiberDim);
_fieldAvg->allocate();
- } else if (_fieldAvg->size() != cells->size()*fiberDim) {
+ } else if (_fieldAvg->chartSize() != cells->size()*fiberDim) {
_fieldAvg->clear();
_fieldAvg->newSection(sectionIn->getChart(), fiberDim);
_fieldAvg->allocate();
Modified: short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/libsrc/meshio/OutputManager.cc 2009-05-28 02:03:25 UTC (rev 15073)
@@ -221,8 +221,19 @@
if (!_fields->hasField(fieldName.c_str())) {
_fields->add(fieldName.c_str(), fieldIn.label());
field_type& fieldOut = _fields->get(fieldName.c_str());
- fieldOut.newSection(fieldIn);
+ const ALE::Obj<typename mesh_type::RealSection>& sectionIn = fieldIn.section();
+ assert(!sectionIn.isNull());
+ const typename mesh_type::RealSection::chart_type chart = sectionIn->getChart();
+ const typename mesh_type::RealSection::chart_type::const_iterator chartBegin =
+ chart.begin();
+ const typename mesh_type::RealSection::chart_type::const_iterator chartEnd =
+ chart.end();
+ const int fiberDim = (chartBegin != chartEnd) ?
+ sectionIn->getFiberDimension(*chartBegin) : 0;
+ fieldOut.newSection(chart, fiberDim);
fieldOut.allocate();
+ fieldOut.vectorFieldType(fieldIn.vectorFieldType());
+ fieldOut.scale(fieldIn.scale());
} // if
field_type& fieldOut = _fields->get(fieldName.c_str());
fieldOut.copy(fieldIn);
@@ -233,8 +244,7 @@
} // if/else
// Satisfy return value. Should never get this far.
- assert(0 != _fields);
- return _fields->get("buffer (dimensioned)");
+ return fieldIn;
} // _dimension
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.cc 2009-05-28 02:03:25 UTC (rev 15073)
@@ -75,22 +75,22 @@
} // spaceDim
// ----------------------------------------------------------------------
-// Get the chart size
+// Get the chart size.
template<typename mesh_type>
int
pylith::topology::Field<mesh_type>::chartSize(void) const
-{ // spaceDim
+{ // chartSize
return _section.isNull() ? 0 : _section->getChart().size();
-} // spaceDim
+} // chartSize
// ----------------------------------------------------------------------
-// Get the size
+// Get the number of degrees of freedom.
template<typename mesh_type>
int
-pylith::topology::Field<mesh_type>::size(void) const
-{ // spaceDim
+pylith::topology::Field<mesh_type>::sectionSize(void) const
+{ // sectionSize
return _section.isNull() ? 0 : _section->size();
-} // spaceDim
+} // sectionSize
// ----------------------------------------------------------------------
// Create seive section.
@@ -211,6 +211,8 @@
std::cout << "Making Field " << _label << " section type 3" << std::endl;
logger.stagePush("Field");
_vecFieldType = src._vecFieldType;
+ _scale = src._scale;
+ _dimensionsOkay = false;
const ALE::Obj<RealSection>& srcSection = src.section();
if (!srcSection.isNull() && _section.isNull()) {
@@ -301,8 +303,8 @@
pylith::topology::Field<mesh_type>::copy(const Field& field)
{ // copy
// Check compatibility of sections
- const int srcSize = (!field._section.isNull()) ? field._section->size() : 0;
- const int dstSize = (!_section.isNull()) ? _section->size() : 0;
+ const int srcSize = field.chartSize();
+ const int dstSize = chartSize();
if (field.spaceDim() != spaceDim() ||
field._vecFieldType != _vecFieldType ||
srcSize != dstSize) {
@@ -351,8 +353,8 @@
pylith::topology::Field<mesh_type>::copy(const ALE::Obj<typename mesh_type::RealSection>& osection)
{ // copy
// Check compatibility of sections
- const int srcSize = (!osection.isNull()) ? osection->size() : 0;
- const int dstSize = (!_section.isNull()) ? _section->size() : 0;
+ const int srcSize = osection->getChart().size();
+ const int dstSize = chartSize();
if (srcSize != dstSize) {
std::ostringstream msg;
@@ -392,8 +394,8 @@
pylith::topology::Field<mesh_type>::operator+=(const Field& field)
{ // operator+=
// Check compatibility of sections
- const int srcSize = (!field._section.isNull()) ? field._section->size() : 0;
- const int dstSize = (!_section.isNull()) ? _section->size() : 0;
+ const int srcSize = field.chartSize();
+ const int dstSize = chartSize();
if (field.spaceDim() != spaceDim() ||
field._vecFieldType != _vecFieldType ||
field._scale != _scale ||
Modified: short/3D/PyLith/trunk/libsrc/topology/Field.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/libsrc/topology/Field.hh 2009-05-28 02:03:25 UTC (rev 15073)
@@ -142,7 +142,7 @@
*
* @returns the number of degrees of freedom.
*/
- int size(void) const;
+ int sectionSize(void) const;
/// Create sieve section.
void newSection(void);
Modified: short/3D/PyLith/trunk/modulesrc/topology/Field.i
===================================================================
--- short/3D/PyLith/trunk/modulesrc/topology/Field.i 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/modulesrc/topology/Field.i 2009-05-28 02:03:25 UTC (rev 15073)
@@ -131,7 +131,7 @@
*
* @returns the number of degrees of freedom.
*/
- int size(void) const;
+ int sectionSize(void) const;
/// Create sieve section.
void newSection(void);
Modified: short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubMesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubMesh.py 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubMesh.py 2009-05-28 02:03:25 UTC (rev 15073)
@@ -46,7 +46,7 @@
ModuleDataWriterVTK.filename(self, self.filename)
ModuleDataWriterVTK.timeFormat(self, self.timeFormat)
- ModuleDataWriterVTK.timeConstant(self, self.timeConstant.value)
+ ModuleDataWriterVTK.timeConstant(self, self.timeConstantN)
return
Modified: short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubSubMesh.py
===================================================================
--- short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubSubMesh.py 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/pylith/meshio/DataWriterVTKSubSubMesh.py 2009-05-28 02:03:25 UTC (rev 15073)
@@ -46,7 +46,7 @@
ModuleDataWriterVTK.filename(self, self.filename)
ModuleDataWriterVTK.timeFormat(self, self.timeFormat)
- ModuleDataWriterVTK.timeConstant(self, self.timeConstant.value)
+ ModuleDataWriterVTK.timeConstant(self, self.timeConstantN)
return
Modified: short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py
===================================================================
--- short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/pylith/perf/MemoryLogger.py 2009-05-28 02:03:25 UTC (rev 15073)
@@ -119,7 +119,7 @@
import pylith.perf.Field
if not stage in self.memory: self.memory[stage] = {}
- fieldModel = pylith.perf.Field.Field(field.label(), field.size(), field.chartSize())
+ fieldModel = pylith.perf.Field.Field(field.label(), field.sectionSize(), field.chartSize())
fieldModel.tabulate(self.memory[stage])
if self.verbose: self.show()
return
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2009-05-28 02:03:25 UTC (rev 15073)
@@ -47,8 +47,6 @@
"""
Formulation.__init__(self, name)
self._loggingPrefix = "TSEx "
- self.solnField = {'name': "disp(t)",
- 'label': "displacement"}
return
@@ -70,30 +68,18 @@
Formulation.initialize(self, dimension, normalizer)
self._info.log("Creating other fields and matrices.")
- self.fields.add("dispIncr(t->t+dt)", "displacement_increment")
self.fields.add("disp(t-dt)", "displacement")
- self.fields.add("residual", "residual")
- self.fields.copyLayout("disp(t)")
- self.fields.solutionName("dispIncr(t->t+dt)")
+ self.fields.copyLayout("dispIncr(t->t+dt)")
self._debug.log(resourceUsageString())
- # Set fields to zero
- lengthScale = normalizer.lengthScale()
- dispIncr = self.fields.get("dispIncr(t->t+dt)")
- dispIncr.scale(lengthScale.value)
- dispIncr.zero()
- dispT = self.fields.get("disp(t)")
- dispT.scale(lengthScale.value)
- dispT.zero()
+ # Setup fields and set to zero
dispTmdt = self.fields.get("disp(t-dt)")
- dispTmdt.scale(lengthScale.value)
dispTmdt.zero()
+ dispT = self.fields.get("disp(t)")
+ dispT.zero()
residual = self.fields.get("residual")
- residual.scale(lengthScale.value)
residual.zero()
- # Create Petsc vectors for fields involved in solve
residual.createVector()
- dispIncr.createVector()
self._debug.log(resourceUsageString())
self._info.log("Creating Jacobian matrix.")
Modified: short/3D/PyLith/trunk/pylith/problems/Formulation.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Formulation.py 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/pylith/problems/Formulation.py 2009-05-28 02:03:25 UTC (rev 15073)
@@ -207,20 +207,35 @@
output.open(totalTime, numTimeSteps)
self._debug.log(resourceUsageString())
+ # Setup fields
self._info.log("Creating solution field.")
- solnName = self.solnField['name']
- self.fields.add(solnName, self.solnField['label'])
- self.fields.solutionName(solnName)
- solution = self.fields.solution()
+ self.fields.add("dispIncr(t->t+dt)", "displacement_increment")
+ self.fields.add("disp(t)", "displacement")
+ self.fields.add("residual", "residual")
+ self.fields.solutionName("dispIncr(t->t+dt)")
+
+ lengthScale = normalizer.lengthScale()
+ solution = self.fields.get("dispIncr(t->t+dt)")
solution.vectorFieldType(solution.VECTOR)
- solution.scale(normalizer.lengthScale().value)
+ solution.scale(lengthScale.value)
solution.newSection(solution.VERTICES_FIELD, dimension)
for constraint in self.constraints:
constraint.setConstraintSizes(solution)
solution.allocate()
for constraint in self.constraints:
constraint.setConstraints(solution)
+ solution = self.fields.get("dispIncr(t->t+dt)")
+ solution.createVector()
solution.createScatter()
+
+ dispT = self.fields.get("disp(t)")
+ dispT.vectorFieldType(dispT.VECTOR)
+ dispT.scale(lengthScale.value)
+
+ residual = self.fields.get("residual")
+ residual.vectorFieldType(residual.VECTOR)
+ residual.scale(lengthScale.value)
+
self._debug.log(resourceUsageString())
self._logger.eventEnd(logEvent)
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2009-05-28 01:50:40 UTC (rev 15072)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2009-05-28 02:03:25 UTC (rev 15073)
@@ -70,8 +70,6 @@
"""
Formulation.__init__(self, name)
self._loggingPrefix = "TSIm "
- self.solnField = {'name': "disp(t)",
- 'label': "displacement"}
self._stepCount = None
return
@@ -90,26 +88,16 @@
"""
Formulation.initialize(self, dimension, normalizer)
+ # Allocate other fields, reusing layout from dispIncr
self._info.log("Creating other fields.")
- self.fields.add("dispIncr(t->t+dt)", "displacement_increment")
- self.fields.add("residual", "residual")
- self.fields.copyLayout("disp(t)")
- self.fields.solutionName("dispIncr(t->t+dt)")
+ self.fields.copyLayout("dispIncr(t->t+dt)")
- # Set fields to zero
- lengthScale = normalizer.lengthScale()
- disp = self.fields.get("disp(t)")
- disp.scale(lengthScale.value)
- disp.zero()
- dispIncr = self.fields.get("dispIncr(t->t+dt)")
- dispIncr.scale(lengthScale.value)
- dispIncr.zero()
+ # Setup fields and set to zero
+ dispT = self.fields.get("disp(t)")
+ dispT.zero()
residual = self.fields.get("residual")
- residual.scale(lengthScale.value)
residual.zero()
- # Create Petsc vectors for fields involved in solve
residual.createVector()
- dispIncr.createVector()
self._debug.log(resourceUsageString())
self._info.log("Creating Jacobian matrix.")
@@ -167,7 +155,6 @@
for constraint in self.constraints:
constraint.setFieldIncr(t, t+dt, dispIncr)
- ### NONLINEAR: Might want to move logic into IntegrateJacobian() and set a flag instead
for integrator in self.integratorsMesh + self.integratorsSubMesh:
integrator.timeStep(dt)
if integrator.needNewJacobian():
@@ -185,22 +172,14 @@
dispIncr = self.fields.get("dispIncr(t->t+dt)")
dispIncr.zero()
- ### NONLINEAR: This moves under SNES control as IntegrateResidual()
- ### NONLINEAR: Also move updateState() from Integrator.poststep() to this function
self._reformResidual(t+dt, dt)
self._info.log("Solving equations.")
residual = self.fields.get("residual")
self._logger.stagePush("Solve")
- #residual.view("RESIDUAL BEFORE SOLVE")
- #self.jacobian.view()
self.solver.solve(dispIncr, self.jacobian, residual)
self._logger.stagePop()
- # BEGIN TEMPORARY
- #dispIncr.view("DISPINCR SOLUTION")
- #residual.view("RESIDUAL")
- # END TEMPORARY
return
@@ -212,6 +191,7 @@
dispIncr = self.fields.get("dispIncr(t->t+dt)")
disp = self.fields.get("disp(t)")
disp += dispIncr
+ dispIncr.zero()
Formulation.poststep(self, t, dt)
More information about the CIG-COMMITS
mailing list