How to compute the force equivalent to the application of a Dirichlet boundary condition

[thelfer]

Dear all,

I have written a small library based on MFEM for non linear mechanics based on a classed derived with NonLinearForm.

I have a small question regarding the postprocessing of my computations.

I impose a displacement on a boundary of my mesh.
I want to compute the force associated with this boundary conditions.

In other codes where boundary conditions are handled by Lagrange multipliers, I just sum up the values of the Lagrange multipliers. I wonder how to perform this computation with MFEM.

Thanks for any help.

[rcarson3]

@thelfer you might find #1842 of interest, since it allows you to apply BCs using constraints and you can get out the lagrange multipliers out from the solver. So, you could get out the forces using the methods you're used to.

@jamiebramwell or @barker29 potentially might have some suggestions as well.

[thelfer]

@rcarson3 Thanks for the link. I currently plan to build the inner forces on the elements adjacent to the boundary and only keep the forces on the boundary (I only need the resultant).
Of course, I am open to any orther ideas !

[thelfer]

I finally managed to compute what I wanted which works whatever the way boundary conditions are handled.

I post here a pseudo-code which may help other users and that can be used in many problems (solid mechanics, heat transfer, etc...). The code that I reproduce here is adapted from the mfem-mgis project, so this code may not compile du to adaptations errors but the overall idea shall be ok.

@rcarson3, @jamiebramwell, @barker29 I would appreciate any comments on this code.

Let:

• m: be a mesh
• fes: be a finite element space
• bid: be a boundary identifier
• i: be a non linear form integrator which computes the inner forces
• u: be the unknowns

Then the resultant force can be computed as follows:

const auto faces = buildFacesDescription(m, fes, bid);
mfem::Vector F;
computeResultantForceOnBoundary(F, fes, i, u, faces);

where :

• buildFacesDescription returns a vector associating to an element of the boundary its adjacent element.
• computeResultantForceOnBoundary actually does the computation

std::vector<std::pair<size_type, size_type>>
buildFacesDescription(mfem::mesh &m, mfem::FiniteElementSpace &fes,
                      const size_type bid) {
  std::vector<std::pair<size_type, size_type>> faces;
  for (size_type i = 0; i < fes.GetNBE(); i++) {
    const auto bdr_attr = m.GetBdrAttribute(i);
    if (bdr_attr != bid) {
      continue;
    }
    const auto *const tr = m.GetBdrFaceTransformations(i);
    if (tr != nullptr) {
      faces.push_back({i, tr->Elem1No});
    }
  }
  return faces;
} // end of buildFaces

void computeResultantForceOnBoundary(
    mfem::Vector &F, mfem::FiniteElementSpace &fes,
    mfem::NonLinearFormIntegrator &i, const mfem::Vector &u,
    const std::vector<std::pair<size_type, size_type>> &faces) {
  const auto nc = fes.GetVDim();
  //   const Vector &px = Prolongate(x); ?? required ??
  mfem::Vector cell_unknowns;
  mfem::Array<int> cell_dofs;
  mfem::Array<int> face_dofs;
  mfem::Vector cell_forces;
  mfem::Vector face_forces;
  F.SetSize(nc);
  F = real{0};
  for (const auto &f : faces) {
    const auto &fe = *(fes.GetFE(std::get<1>(f)));
    auto &tr = *(fes.GetElementTransformation(std::get<1>(f)));
    fes.GetElementVDofs(std::get<1>(f), cell_dofs);
    y.GetSubVector(vdofs, cell_unknowns);
    // compute the inner forces
    i.AssembleElementVector(cell_forces, fe, tr, cell_unknowns);
    // computation of the resultant by summing the contributions on the
    // considered boundary
    fes.GetBdrElementVDofs(std::get<0>(f), face_dofs);
    face_forces.SetSize(face_dofs.Size());
    const auto n_cell_nodes = fe.GetDof();
    const auto n_face_nodes = face_dofs.Size() / fe.GetDim();
    const auto cell_dofs_begin = cell_dofs.GetData();
    // here we take into account the fact that the components are stored
    // contiguously
    for (size_type c = 0; c != nc; ++c) {
      const auto cell_dofs_c_begin = cell_dofs_begin + c * n_cell_nodes;
      const auto cell_dofs_c_end = cell_dofs_c_begin + n_cell_nodes;
      for (size_type i = 0; i != n_face_nodes; ++i) {
        const auto fid = i + c * n_face_nodes;
        const auto face_dof = face_dofs[fid];
        const auto pc = std::find(cell_dofs_c_begin, cell_dofs_c_end, face_dof);
        F[c] += cell_forces[pc - cell_dofs_begin];
      }
    }
  }
} // end of computeResultantForceOnBoundary