Fix boundary integration overlap, parallel output of refined mesh, and a few other small fixes

configure.ac

@@ -114,7 +114,7 @@ AC_CHECK_HEADER([mpi.h], [], [AC_MSG_ERROR([header 'mpi.h' not found])])
 
 dnl PETSC
 AC_LANG(C)
-CIT_PATH_PETSC([3.23.2])
+CIT_PATH_PETSC([3.23.4])
 CIT_HEADER_PETSC
 CIT_CHECK_LIB_PETSC
 

docs/user/components/materials/Elasticity.md

@@ -31,7 +31,7 @@ Implements `Material`.
 * `description`=\<str\>: Descriptive label for material (deprecated in v5.0; will be removed in v6.0).
   - **default value**: ''
   - **current value**: '', from {default}
-  - **validator**: <function validateDescription at 0x10608f2e0>
+  - **validator**: <function validateDescription at 0x10480b2e0>
 * `label`=\<str\>: Name of label for material. Currently only 'material-id' is allowed.
   - **default value**: 'material-id'
   - **current value**: 'material-id', from {default}
@@ -50,7 +50,6 @@ Example of setting `Elasticity` Pyre properties and facilities in a parameter fi
 :::{code-block} cfg
 [pylithapp.problem.materials.mat_elastic]
 label_value = 4
-description = Upper crust elastic material
 use_body_force = False
 bulk_rheology = pylith.materials.IsotropicLinearElasticity
 

docs/user/components/materials/IncompressibleElasticity.md

@@ -31,7 +31,7 @@ Implements `Material`.
 * `description`=\<str\>: Descriptive label for material (deprecated in v5.0; will be removed in v6.0).
   - **default value**: ''
   - **current value**: '', from {default}
-  - **validator**: <function validateDescription at 0x10608f2e0>
+  - **validator**: <function validateDescription at 0x10480b2e0>
 * `label`=\<str\>: Name of label for material. Currently only 'material-id' is allowed.
   - **default value**: 'material-id'
   - **current value**: 'material-id', from {default}
@@ -49,7 +49,6 @@ Example of setting `IncompressibleElasticity` Pyre properties and facilities in
 
 :::{code-block} cfg
 [pylithapp.problem.materials.mat_incompelastic]
-description = Upper crust incompressible elastic material
 label_value = 3
 use_body_force = True
 bulk_rheology = pylith.materials.IsotropicLinearIncompElasticity

docs/user/components/materials/Poroelasticity.md

@@ -31,7 +31,7 @@ Implements `Material`.
 * `description`=\<str\>: Descriptive label for material (deprecated in v5.0; will be removed in v6.0).
   - **default value**: ''
   - **current value**: '', from {default}
-  - **validator**: <function validateDescription at 0x10608f2e0>
+  - **validator**: <function validateDescription at 0x10480b2e0>
 * `label`=\<str\>: Name of label for material. Currently only 'material-id' is allowed.
   - **default value**: 'material-id'
   - **current value**: 'material-id', from {default}
@@ -55,7 +55,6 @@ Example of setting `Poroelasticity` Pyre properties and facilities in a paramete
 
 :::{code-block} cfg
 [pylithapp.problem.materials.mat_poroelastic]
-description = Upper crust poroelastic material
 label_value = 3
 use_body_force = True
 use_source_density = False

docs/user/components/meshio/MeshIOPetsc.md

@@ -24,7 +24,7 @@ Implements `MeshIOObj`.
 * `filename`=\<str\>: Name of mesh file for reading with PETSc.
   - **default value**: ''
   - **current value**: '', from {default}
-* `gmsh_mark_vertices`=\<bool\>: Gmsh file marks faces, edges, and vertices rather than just faces.
+* `gmsh_mark_recursive`=\<bool\>: Gmsh file marks faces, edges, and vertices rather than just faces (3D) or edges (2D).
   - **default value**: False
   - **current value**: False, from {default}
 * `options_prefix`=\<str\>: Name of PETSc options prefix for this mesh.

docs/user/components/meshio/OutputPhysics.md

@@ -45,8 +45,8 @@ Example of setting `OutputPhysics` Pyre properties and facilities in a parameter
 :::{code-block} cfg
 [observer]
 # Skip two time steps between output.
-output_trigger = pylith.meshio.OutputTriggerStep
-output_trigger.num_skip = 2
+trigger = pylith.meshio.OutputTriggerStep
+trigger.num_skip = 2
 
 # Write output to HDF5 file with name `boundary_xpos.h5`.
 writer = pylith.meshio.DataWriterHDF5

docs/user/components/meshio/OutputSolnBoundary.md

@@ -29,7 +29,7 @@ Implements `OutputSoln`.
 * `label`=\<str\>: Name of label identifier for external boundary.
   - **default value**: ''
   - **current value**: '', from {default}
-  - **validator**: <function validateLabel at 0x10816d800>
+  - **validator**: <function validateLabel at 0x104835e40>
 * `label_value`=\<int\>: Value of label identifier for external boundary (tag of physical group in Gmsh files).
   - **default value**: 1
   - **current value**: 1, from {default}
@@ -53,8 +53,8 @@ data_fields = [displacement]
 label = boundary_xpos
 
 # Skip two time steps between output.
-output_trigger = pylith.meshio.OutputTriggerStep
-output_trigger.num_skip = 2
+trigger = pylith.meshio.OutputTriggerStep
+trigger.num_skip = 2
 
 # Write output to HDF5 file with name `boundary_xpos.h5`.
 writer = pylith.meshio.DataWriterHDF5

docs/user/components/meshio/OutputSolnDomain.md

@@ -44,8 +44,8 @@ Example of setting `OutputSolnDomain` Pyre properties and facilities in a parame
 data_fields = [displacement]
 
 # Skip two time steps between output.
-output_trigger = pylith.meshio.OutputTriggerStep
-output_trigger.num_skip = 2
+trigger = pylith.meshio.OutputTriggerStep
+trigger.num_skip = 2
 
 # Write output to HDF5 file with name `domain.h5`.
 writer = pylith.meshio.DataWriterHDF5

docs/user/components/meshio/OutputSolnPoints.md

@@ -55,8 +55,8 @@ reader = pylith.meshio.PointsList
 reader.filename = stations.txt
 
 # Skip two time steps between output.
-output_trigger = pylith.meshio.OutputTriggerStep
-output_trigger.num_skip = 2
+trigger = pylith.meshio.OutputTriggerStep
+trigger.num_skip = 2
 
 # Write output to HDF5 file with name `domain.h5`.
 writer = pylith.meshio.DataWriterHDF5

docs/user/components/topology/MeshImporter.md

@@ -13,7 +13,7 @@ Implements `MeshGenerator`.
 * `distributor`: Distributes mesh among processes.
   - **current value**: 'mesh_distributor', from {default}
   - **configurable as**: mesh_distributor, distributor
-* `reader`: Reader for mesh files.
+* `reader`: Reader for mesh file.
   - **current value**: 'meshioascii', from {default}
   - **configurable as**: meshioascii, reader
 * `refiner`: Performs uniform global mesh refinement after distribution among processes (default is no refinement).

libsrc/pylith/bc/AbsorbingDampers.cc

@@ -121,6 +121,7 @@ pylith::bc::AbsorbingDampers::createIntegrator(const pylith::topology::Field& so
     integrator->setSubfieldName(getSubfieldName());
     integrator->setLabelName(getLabelName());
     integrator->setLabelValue(getLabelValue());
+    integrator->createLabelDS(solution, solution.getMesh().getDimension()-1);
 
     _AbsorbingDampers::setKernelsResidual(integrator, *this, solution);
     BoundaryCondition::_setKernelsDiagnosticField(integrator, solution);
@@ -201,8 +202,8 @@ pylith::bc::_AbsorbingDampers::setKernelsResidual(pylith::feassemble::Integrator
           << "_AbsorbingDampers::_setKernelsRHSResidual(integrator="<<integrator<<", bc="<<typeid(bc).name()
           <<", solution="<<solution.getLabel()<<")"<<pythia::journal::endl;
 
-    PetscBdPointFunc g0 = pylith::fekernels::AbsorbingDampers::g0;
-    PetscBdPointFunc g1 = NULL;
+    PetscBdPointFn* g0 = pylith::fekernels::AbsorbingDampers::g0;
+    PetscBdPointFn* g1 = NULL;
 
     std::vector<ResidualKernels> kernels(1);
     kernels[0] = ResidualKernels(bc.getSubfieldName(), pylith::feassemble::Integrator::RHS, g0, g1);

libsrc/pylith/bc/DirichletTimeDependent.cc

@@ -338,7 +338,7 @@ pylith::bc::_DirichletTimeDependent::setKernelConstraint(pylith::feassemble::Con
           << "setKernelConstraint(constraint="<<constraint<<", bc="<<typeid(bc).name()<<", solution="<<solution.getLabel()
           <<")" << pythia::journal::endl;
 
-    PetscBdPointFunc bcKernel = NULL;
+    PetscBdPointFn* bcKernel = NULL;
 
     const pylith::topology::Field::VectorFieldEnum fieldType = solution.getSubfieldInfo(bc.getSubfieldName()).description.vectorFieldType;
     const bool isScalarField = fieldType == pylith::topology::Field::SCALAR;

libsrc/pylith/bc/NeumannTimeDependent.cc

@@ -201,6 +201,7 @@ pylith::bc::NeumannTimeDependent::createIntegrator(const pylith::topology::Field
     integrator->setSubfieldName(getSubfieldName());
     integrator->setLabelName(getLabelName());
     integrator->setLabelValue(getLabelValue());
+    integrator->createLabelDS(solution, solution.getMesh().getDimension()-1);
 
     _NeumannTimeDependent::setKernelsResidual(integrator, *this, solution, _formulation);
     BoundaryCondition::_setKernelsDiagnosticField(integrator, solution);
@@ -338,8 +339,8 @@ pylith::bc::_NeumannTimeDependent::setKernelsResidual(pylith::feassemble::Integr
     const int bitTimeHistory = bc.useTimeHistory() ? 0x4 : 0x0;
     const int bitUse = bitInitial | bitRate | bitTimeHistory;
 
-    PetscBdPointFunc r0 = NULL;
-    PetscBdPointFunc r1 = NULL;
+    PetscBdPointFn* r0 = NULL;
+    PetscBdPointFn* r1 = NULL;
     switch (bitUse) {
     case 0x1:
         r0 = (isScalarField) ? pylith::fekernels::NeumannTimeDependent::f0_initial_scalar :

libsrc/pylith/bc/NeumannUserFn.cc

@@ -90,7 +90,7 @@ pylith::bc::NeumannUserFn::deallocate(void) {
 // ------------------------------------------------------------------------------------------------
 // Set user function specifying field on boundary.
 void
-pylith::bc::NeumannUserFn::setUserFn(PetscBdPointFunc fn) {
+pylith::bc::NeumannUserFn::setUserFn(PetscBdPointFn* fn) {
     PYLITH_COMPONENT_DEBUG("setUserFn(fn="<<fn<<")");
 
     _fn = fn;
@@ -99,7 +99,7 @@ pylith::bc::NeumannUserFn::setUserFn(PetscBdPointFunc fn) {
 
 // ------------------------------------------------------------------------------------------------
 // Get user function specifying field on boundary.
-PetscBdPointFunc
+PetscBdPointFn*
 pylith::bc::NeumannUserFn::getUserFn(void) const {
     return _fn;
 } // getUserFn
@@ -116,6 +116,7 @@ pylith::bc::NeumannUserFn::createIntegrator(const pylith::topology::Field& solut
     integrator->setSubfieldName(getSubfieldName());
     integrator->setLabelName(getLabelName());
     integrator->setLabelValue(getLabelValue());
+    integrator->createLabelDS(solution, solution.getMesh().getDimension()-1);
 
     _NeumannUserFn::setKernelsResidual(integrator, *this, solution, _formulation);
 
@@ -178,8 +179,8 @@ pylith::bc::_NeumannUserFn::setKernelsResidual(pylith::feassemble::IntegratorBou
           << solution.getLabel()<<")"
           << pythia::journal::endl;
 
-    PetscBdPointFunc r0 = bc.getUserFn();
-    PetscBdPointFunc r1 = NULL;
+    PetscBdPointFn* r0 = bc.getUserFn();
+    PetscBdPointFn* r1 = NULL;
 
     std::vector<ResidualKernels> kernels(1);
     switch (formulation) {

libsrc/pylith/bc/NeumannUserFn.hh

@@ -32,13 +32,13 @@ public:
      *
      * @param[in] fn Function specifying field on boundary.
      */
-    void setUserFn(PetscBdPointFunc fn);
+    void setUserFn(PetscBdPointFn* fn);
 
     /** Get user function specifying field on boundary
      *
      * @preturns Function specifying field on boundary.
      */
-    PetscBdPointFunc getUserFn(void) const;
+    PetscBdPointFn* getUserFn(void) const;
 
     /** Create integrator and set kernels.
      *
@@ -76,7 +76,7 @@ protected:
     // PRIVATE MEMBERS ////////////////////////////////////////////////////////////////////////////
 private:
 
-    PetscBdPointFunc _fn; ///< Kernel for boundary integration.
+    PetscBdPointFn* _fn; ///< Kernel for boundary integration.
 
     // NOT IMPLEMENTED ////////////////////////////////////////////////////////////////////////////
 private:

libsrc/pylith/faults/FaultCohesive.cc

@@ -395,6 +395,7 @@ pylith::faults::FaultCohesive::createIntegrator(const pylith::topology::Field& s
     integrator->setLabelName(getCohesiveLabelName());
     integrator->setLabelValue(getCohesiveLabelValue());
     integrator->setSurfaceLabelName(getSurfaceLabelName());
+    integrator->createLabelDS(solution, solution.getMesh().getDimension()-1);
 
     pylith::feassemble::InterfacePatches* patches =
         pylith::feassemble::InterfacePatches::createMaterialPairs(this, solution.getDM());

libsrc/pylith/faults/FaultCohesiveImpulses.cc

@@ -309,16 +309,16 @@ pylith::faults::FaultCohesiveImpulses::_setKernelsResidual(pylith::feassemble::I
     switch (_formulation) {
     case pylith::problems::Physics::QUASISTATIC: {
         // Elasticity equation (displacement) for negative side of the fault.
-        const PetscBdPointFunc f0u_neg = pylith::fekernels::FaultCohesiveKin::f0u_neg;
-        const PetscBdPointFunc f1u_neg = NULL;
+        PetscBdPointFn* f0u_neg = pylith::fekernels::FaultCohesiveKin::f0u_neg;
+        PetscBdPointFn* f1u_neg = NULL;
 
         // Elasticity equation (displacement) for positive side of the fault.
-        const PetscBdPointFunc f0u_pos = pylith::fekernels::FaultCohesiveKin::f0u_pos;
-        const PetscBdPointFunc f1u_pos = NULL;
+        PetscBdPointFn* f0u_pos = pylith::fekernels::FaultCohesiveKin::f0u_pos;
+        PetscBdPointFn* f1u_pos = NULL;
 
         // Fault slip constraint equation.
-        const PetscBdPointFunc f0l = pylith::fekernels::FaultCohesiveKin::f0l_slip;
-        const PetscBdPointFunc f1l = NULL;
+        PetscBdPointFn* f0l = pylith::fekernels::FaultCohesiveKin::f0l_slip;
+        PetscBdPointFn* f1l = NULL;
 
         kernels.resize(3);
         kernels[0] = ResidualKernels("displacement", integrator_t::LHS, integrator_t::NEGATIVE_FACE,
@@ -357,20 +357,20 @@ pylith::faults::FaultCohesiveImpulses::_setKernelsJacobian(pylith::feassemble::I
     std::vector<JacobianKernels> kernels;
     switch (_formulation) {
     case QUASISTATIC: {
-        const PetscBdPointJac Jf0ul_neg = pylith::fekernels::FaultCohesiveKin::Jf0ul_neg;
-        const PetscBdPointJac Jf1ul_neg = NULL;
-        const PetscBdPointJac Jf2ul_neg = NULL;
-        const PetscBdPointJac Jf3ul_neg = NULL;
-
-        const PetscBdPointJac Jf0ul_pos = pylith::fekernels::FaultCohesiveKin::Jf0ul_pos;
-        const PetscBdPointJac Jf1ul_pos = NULL;
-        const PetscBdPointJac Jf2ul_pos = NULL;
-        const PetscBdPointJac Jf3ul_pos = NULL;
-
-        const PetscBdPointJac Jf0lu = pylith::fekernels::FaultCohesiveKin::Jf0lu;
-        const PetscBdPointJac Jf1lu = NULL;
-        const PetscBdPointJac Jf2lu = NULL;
-        const PetscBdPointJac Jf3lu = NULL;
+        PetscBdPointJacFn* Jf0ul_neg = pylith::fekernels::FaultCohesiveKin::Jf0ul_neg;
+        PetscBdPointJacFn* Jf1ul_neg = NULL;
+        PetscBdPointJacFn* Jf2ul_neg = NULL;
+        PetscBdPointJacFn* Jf3ul_neg = NULL;
+
+        PetscBdPointJacFn* Jf0ul_pos = pylith::fekernels::FaultCohesiveKin::Jf0ul_pos;
+        PetscBdPointJacFn* Jf1ul_pos = NULL;
+        PetscBdPointJacFn* Jf2ul_pos = NULL;
+        PetscBdPointJacFn* Jf3ul_pos = NULL;
+
+        PetscBdPointJacFn* Jf0lu = pylith::fekernels::FaultCohesiveKin::Jf0lu;
+        PetscBdPointJacFn* Jf1lu = NULL;
+        PetscBdPointJacFn* Jf2lu = NULL;
+        PetscBdPointJacFn* Jf3lu = NULL;
 
         kernels.resize(3);
         const char* nameDisplacement = "displacement";