[Boundary Condition] Nonconforming traction BC (copy from xmpm branch) by jgiven100 · Pull Request #66 · geomechanics/mpm

jgiven100 · 2022-08-17T21:09:53Z

Describe the PR

Overview

Add new class for NonconformingTractionConstraint that allows for a Nonconforming Tractiong Constraint to be applied. The constraint supports either constant pressure (e.g., pressure vessel problem) or hydrostatic pressure (e.g., levee probem) along a free surface. These features are toggled via input.JSON.

For nodes "close" to the boundary, f_int and f_ext must be updated.

Internal nodal force is updated such that:

Add applied stress at particle level
Subtract applied stress at bounary at cell level

External nodal force is updated such that:

Add divergence of stress at particle level
Subtract divergence of stress at cell level

These addition/subtraction rules are based on imposing a virtual stress field that produces an exactly equivalent resposne within the material to a traditional boundary traction.

Simple Boundary Detection

PR Includes a simple algorithm to detect which cells are along the material boundary.

    // Find cells and nodes located at the vicinity of the interface
    for (auto citr = cells_.cbegin(); citr != cells_.cend(); citr++) {
      // Loop over void cells
      if ((*citr)->nparticles() != 0) continue;

      // Find non-void cell connected with void cell
      for (auto cell_neigh : (*citr)->neighbours()) {
        if (map_cells_[cell_neigh]->nparticles() == 0) continue;
        boundary_cell_list.insert((*citr)->id());
        boundary_cell_list.insert(map_cells_[cell_neigh]->id());
        for (auto node : map_cells_[cell_neigh]->nodes())
          boundary_node_list.insert(node->id());
      }
    }

Therefore, only nodes located "close" to the boundary (and have particles within their support) have f_int and f_ext updated.

Related Issues/PRs
N/a.

Additional context
N/a.

jgiven100 · 2022-08-20T17:08:51Z

Currently there are 2 options

Overall notes

More than 1 "nonconforming_pressure_constraints" can be defined.
Optional "inside" boolean flag for bounding box to indicate if the non-conforming surface is inside or outside the bounding box. Default is true.

Constant Pressure

For the constant pressure scenario (e.g. thick-walled cylinder), the following inputs are required:

  "mesh": {
    ...,
    "boundary_conditions": {
      ...,
      "nonconforming_traction_constraints": [
        {
          "bounding_box": [-3, 3, -3, 3, -1, 1],
          "inside": true,
          "math_function_id": 0,
          "pressure": 1e+6
        }
      ]
    }

Specific Notes:

"math_function_id" is optional. If no id is provided, the pressure will be instantaneously applied.

Hydrostatic Pressure

For the hydrostatic pressure scearnio (e.g. body in fluid), the following inputs are required:

  "mesh": {
    ...,
    "boundary_conditions": {
      ...,
      "nonconforming_traction_constraints": [
        {
          "bounding_box": [-3, 3, -3, 3, -1, 1],
          "datum": 1.0,
          "fluid_density": 1000,
          "hydrostatic": true
        }
      ]
    }, 
  "external_loading_conditions": {
    "gravity": [0, -10]
  },

Specific Notes

"gravity" must be provided in "external_loading_conditions".
Hydrostatic pressure acts in y-dir for 2d case and z-dir for 3d case.

codecov · 2022-08-20T22:33:32Z

Codecov Report

Attention: Patch coverage is 75.89744% with 94 lines in your changes missing coverage. Please review.

Project coverage is 95.56%. Comparing base (991be81) to head (1249585).
Report is 15 commits behind head on master.

Files with missing lines	Patch %	Lines
include/solvers/mpm_base.tcc	20.00%	44 Missing ⚠️
...lude/loads_bcs/nonconforming_traction_constraint.h	43.86%	32 Missing ⚠️
include/mesh/mesh.tcc	91.43%	9 Missing ⚠️
include/cells/cell.tcc	76.92%	6 Missing ⚠️
include/particles/particle_base.h	0.00%	2 Missing ⚠️
include/solvers/mpm_explicit.tcc	50.00%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master      #66      +/-   ##
==========================================
- Coverage   95.71%   95.56%   -0.15%     
==========================================
  Files         241      242       +1     
  Lines       50259    50650     +391     
==========================================
+ Hits        48102    48399     +297     
- Misses       2157     2251      +94

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

yliang-sn

It looks good to me. Thanks!

include/loads_bcs/nonconforming_pressure_constraint.h

include/particles/particle.h

bodhinandach · 2022-08-29T21:03:40Z

@jgiven100 Is this PR ready for review?

jgiven100 · 2022-08-29T21:05:21Z

@jgiven100 Is this PR ready for review?

Yes.

bodhinandach · 2022-08-29T21:11:23Z

Can you change the PR description then to help me understand the design?

jgiven100 · 2022-09-01T00:07:36Z

Will pull all updates from master at 1 time after #65 is merged

bodhinandach · 2022-09-01T00:48:15Z

I think you can pull now. I want to check #65 again later.

bodhinandach

@yliang-sn @jgiven100 Thank you very much for your effort. This looks great. I have some comments which can improve the code. Also, I can help to make this parallel. The only bottleneck currently is to apply the forces from mesh, which I noted in my comments. Let's discuss this in person.

bodhinandach · 2022-09-17T01:13:24Z

include/loads_bcs/nonconforming_traction_constraint.h

+  // Fluid density for hydrostatic free surface
+  double fluid_density() const { return fluid_density_; }


I think this should be called material_density instead of fluid density.

bodhinandach · 2022-09-17T01:14:50Z

include/mesh/mesh.h

+  bool create_nonconforming_traction_constraint(
+      const std::vector<double> bounding_box, const double datum,
+      const double fluid_density, const double gravity, const bool hydrostatic,
+      const bool inside, const std::shared_ptr<FunctionBase>& mfunction,
+      const double pressure);


Can we use "assign" instead of "create" to unify it with other constraints?

bodhinandach · 2022-09-17T01:32:41Z

include/mesh/mesh.tcc

+      // Find non-void cell connected with void cell
+      for (auto cell_neigh : (*citr)->neighbours()) {
+        if (map_cells_[cell_neigh]->nparticles() == 0) continue;
+        boundary_cell_list.insert((*citr)->id());
+        boundary_cell_list.insert(map_cells_[cell_neigh]->id());
+        for (auto node : map_cells_[cell_neigh]->nodes())
+          boundary_node_list.insert(node->id());
+      }


This won't work in MPI

bodhinandach · 2022-09-17T01:33:55Z

include/mesh/mesh.tcc

+    Eigen::Matrix<double, 6, 1> traction;
+    traction.setZero();


This should be stress right? Since it is a 2nd order tensor. Traction is a vector.

bodhinandach · 2022-09-17T01:35:18Z

include/mesh/mesh.tcc

+      // // Placeholder for arbitrary traction and divergence terms
+      // else {
+      //   traction << 0., 0., 0., 0., 0., 0.;
+      //   divergence_traction << 0., 0., 0.;
+      // }


How about defining this as an option in the initial input file?

bodhinandach · 2022-09-17T01:50:31Z

include/particles/particle.h

+  //! Minus the internal force of the virtual stress field
+  //! \param[in] traction Boundary traction
+  //! \param[in] divergence_traction Divergence of boundary traction
+  virtual void minus_virtual_stress_field(


update or apply instead of minus?

bodhinandach · 2022-09-17T01:56:23Z

include/particles/particle.tcc

+    if (Tdim == 2) {
+      force[0] = dn_dx_(i, 0) * traction[0] + dn_dx_(i, 1) * traction[3];
+      force[1] = dn_dx_(i, 0) * traction[3] + dn_dx_(i, 1) * traction[1];
+    } else if (Tdim == 3) {
+      force[0] = dn_dx_(i, 0) * traction[0] + dn_dx_(i, 1) * traction[3] +
+                 dn_dx_(i, 2) * traction[5];
+      force[1] = dn_dx_(i, 0) * traction[3] + dn_dx_(i, 1) * traction[1] +
+                 dn_dx_(i, 2) * traction[4];
+      force[2] = dn_dx_(i, 0) * traction[5] + dn_dx_(i, 1) * traction[4] +
+                 dn_dx_(i, 2) * traction[2];
+    }


I would suggest a template specialization so that you dont have to check the if else condition every time step.

bodhinandach · 2022-09-17T01:59:15Z

include/particles/particle.tcc

+    for (unsigned j = 0; j < Tdim; j++)
+      force[j] += shapefn_[i] * divergence_traction[j] * this->volume_;


Why can't you do it without the loop?

Suggested change

for (unsigned j = 0; j < Tdim; j++)

force[j] += shapefn_[i] * divergence_traction[j] * this->volume_;

force.noalias() += shapefn_[i] * divergence_traction * this->volume_;

jgiven100 · 2024-10-01T20:36:56Z

@bodhinandach Branch is successfully rebased:

All tests passing
I can reproduce previous results (locally)

Next steps:

Update algorithm for nonlocal
Increase coverage

bodhinandach · 2024-10-01T22:56:29Z

Thanks @jgiven100 I will have a look into this. Does this utilize the optimized algorithm that you wrote in CPM or does this still use the initial approach written in CMAME?

jgiven100 · 2024-10-01T23:11:23Z

This corresponds to the initial approach as described in CMAME paper.

If we use the "extra efficient" approach described in the CPM paper, I would not recommend coupling with high-order shape function. The non-constant shape function gradients have too much sensitivity on particle position, which introduces unacceptable instability.

Since you have mentioned you are most interested in using VSB method + B-Spline MPM, we can very easily extend the CMAME approach.

bodhinandach · 2024-10-01T23:35:18Z

Understood, thanks!

jgiven100 requested a review from yliang-sn August 17, 2022 21:09

jgiven100 self-assigned this Aug 17, 2022

This comment was marked as outdated.

Sign in to view

jgiven100 added the enhancement New feature or request label Aug 23, 2022

jgiven100 requested a review from bodhinandach August 23, 2022 22:26

jgiven100 marked this pull request as ready for review August 23, 2022 22:27

jgiven100 changed the title ~~copy nonconforming pressure BC from xmpm branch~~ [Boundary Condition] Nonconforming pressure BC (copy from xmpm branch) Aug 23, 2022

yliang-sn previously approved these changes Aug 29, 2022

View reviewed changes

include/loads_bcs/nonconforming_pressure_constraint.h Outdated Show resolved Hide resolved

include/particles/particle.h Show resolved Hide resolved

jgiven100 changed the title ~~[Boundary Condition] Nonconforming pressure BC (copy from xmpm branch)~~ [Boundary Condition] Nonconforming traction BC (copy from xmpm branch) Aug 29, 2022

jgiven100 dismissed yliang-sn’s stale review via fd701af August 29, 2022 22:20

bodhinandach requested changes Sep 17, 2022

View reviewed changes

bodhinandach marked this pull request as draft December 6, 2022 00:12

jgiven100 and others added 10 commits October 1, 2024 10:49

copy nonconforming pressure BC from xmpm branch

80910b9

hydrostatic add, not working

97b53c9

add support for hydrostatic pressure

25d285c

fix 2d, 3d, Tdim issues

a63fa3a

shift constant pressure in code

debf49f

start to add testing

b2a9b59

increase coverage

d06fadb

increase coverage pt 2

8f13890

the great renaming... pressure -> traction

f8a701d

compute the particle pressure for hydrostatic conditions

c7b70b3

jgiven100 force-pushed the non-conforming-pressure branch from 3451dae to c7b70b3 Compare October 1, 2024 17:50

tidy after rebase

bd99402

make vsb method compatible with non-local shape function

1249585

jgiven100 force-pushed the non-conforming-pressure branch from 860eec6 to 1249585 Compare October 10, 2024 16:24

		// Fluid density for hydrostatic free surface
		double fluid_density() const { return fluid_density_; }

		for (unsigned j = 0; j < Tdim; j++)
		force[j] += shapefn_[i] * divergence_traction[j] * this->volume_;

	for (unsigned j = 0; j < Tdim; j++)
	force[j] += shapefn_[i] * divergence_traction[j] * this->volume_;
	force.noalias() += shapefn_[i] * divergence_traction * this->volume_;

Conversation

jgiven100 commented Aug 17, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Overview

Simple Boundary Detection

Uh oh!

This comment was marked as outdated.

This comment was marked as outdated.

jgiven100 commented Aug 20, 2022 • edited by LDTalbot Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Currently there are 2 options

Overall notes

Constant Pressure

Specific Notes:

Hydrostatic Pressure

Specific Notes

Uh oh!

codecov bot commented Aug 20, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

yliang-sn left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

bodhinandach commented Aug 29, 2022

Uh oh!

jgiven100 commented Aug 29, 2022

Uh oh!

bodhinandach commented Aug 29, 2022

Uh oh!

jgiven100 commented Sep 1, 2022

Uh oh!

bodhinandach commented Sep 1, 2022

Uh oh!

bodhinandach left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

jgiven100 commented Oct 1, 2024

Uh oh!

bodhinandach commented Oct 1, 2024

Uh oh!

jgiven100 commented Oct 1, 2024

Uh oh!

bodhinandach commented Oct 1, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

jgiven100 commented Aug 17, 2022 •

edited

Loading

jgiven100 commented Aug 20, 2022 •

edited by LDTalbot

Loading

codecov bot commented Aug 20, 2022 •

edited

Loading