Projecting 2D mechanical fields onto 1D lower-dimensional elements

[Jensoulou]

Hi,

I am working with lower-dimensional elements (1D embedded in a 2D domain). I want to apply a tensile deformation perpendicular to the direction of the fracture modeled by the lower-dimensional elements (see Figure). The fracture should be mechanically invisible, i.e. the two sides of the fracture remain bonded and do not open or slip.

My goal is to use the 2D mechanical information (displacements, strains, stresses, etc.) within the lower-dimensional (1D) elements.

Question: How can I project or otherwise access the bulk mechanical fields from the 2D elements onto the 1D lower-dimensional elements in MOOSE?

Many thanks,

Jens

[GiudGiud]

@lindsayad ?

If all the 1D elements are grouped into a sideset, objects such as integrated boundary conditions can be executed there with the interior parent (the 2D element) material properties being accessible.

Alternatively, you can look at the mortar constraint code, which also relies on these lowerD elements

[lindsayad]

What do you want to do on the 1D elements? Do you have variables block-restricted there and want to fill residuals for them? With this conforming mesh, some kind of interface kernel might be the best although it will certainly need to override computeResidual() and computeJacobian() since you are not interested in filling residuals for the higher-dimensional element degrees of freedom. But like mortar, an interface kernel could give you access to higher-d elem and neighbor side variables values and material computations as well as access to variable values from the lower-d elements.

This would need a change in the framework paths for interface kernel computations (right now lower-D elements are not reinit'd on that path) but the fix would not be hard

[Jensoulou]

Thank you very much for both your answers!
Here is how I currently understand things - please correct me if I’m misunderstanding:

• Mortar code:
  I’ve looked at the mortar constraints classes. They seem designed to impose values from one side subdomain to the other side subdomain through the mortar interface. I don’t really see a way there to compute values inside the mortar elements themselves.

• Interface kernels:
  Interface kernels appear to be built on the neighbor material property framework. But lower-dimensional elements don’t have higher-dimensional elements as neighbors, since “neighbors” in MOOSE share a side. A side of a lower-dimensional element is two dimensions smaller than the bulk higher-D element, so they cannot directly share a side. This makes the higher-D element invisible to the lower-D element in the neighbor sense.

• Hijacking the path:
  Actually, it is not totally invisible to each other because there is the 'interior_parent' information. But there is no material data framework using the interior parent, only these options

enum MaterialDataType
{
  BLOCK_MATERIAL_DATA,
  BOUNDARY_MATERIAL_DATA,
  FACE_MATERIAL_DATA,
  NEIGHBOR_MATERIAL_DATA,
  INTERFACE_MATERIAL_DATA
};

Even if lower-D elements are reinitialized along the interface path, there is no way to retrieve the necessary material properties from the higher-D parent element.

My concrete goal is to add mechanical couplings in a porous flow model with faults as lower-dimensional elements. For example, in PorousFlowMassVolumetricExpansion, the strain rate is used. In a lower-D element, the computation of strain rate loses the contribution of the missing direction, so the volumetric expansion does not match the bulk computation at the same location. To be consistent, Ι believe the fracture elements need to either (a) recover the missing directional information from the bulk, or (b) directly project the bulk strain rate (or stress, displacement, etc.) onto their quadrature points.

Is my goal and my problem clearer to you? Do you have any suggestions that might help me solve it?

[lindsayad]

It's still not clear to me what you want to compute on the lower-d elements sorry. I can be a little dense. Is it material properties that you would then use in the higher-D element computations? It doesn't sound like you are adding any variables restricted to the lower-D elements?

We definitely use higer-D face computations within mortar constraints which are "operating" on lower-dimensional mortar segments and we allow pulling in lower-d variables into higher-D face material computations. Perhaps this code excerpt from loopOverMortarSegments in MortarUtils.h would be helpful to see

      // reinit the variables/residuals/jacobians on the secondary interior
      fe_problem.reinitElemFaceRef(
          reinit_secondary_elem, secondary_side_id, TOLERANCE, &xi1_pts, nullptr, tid);

      const Elem * reinit_primary_elem = primary_ip;

      // If we're on the displaced mesh, we need to get the corresponding undisplaced elem before
      // calling fe_problem.reinitElemFaceRef
      if (displaced)
        reinit_primary_elem = fe_problem.mesh().elemPtr(reinit_primary_elem->id());

      // reinit the variables/residuals/jacobians on the primary interior
      fe_problem.reinitNeighborFaceRef(
          reinit_primary_elem, primary_side_id, TOLERANCE, &xi2_pts, nullptr, tid);

      // reinit neighbor materials, but be careful not to execute stateful materials since
      // conceptually they don't make sense with mortar (they're not interpolary)
      fe_problem.reinitMaterialsNeighbor(primary_ip->subdomain_id(),
                                         /*tid=*/tid,
                                         /*swap_stateful=*/false,
                                         &primary_ip_mats);

      // reinit the variables/residuals/jacobians on the lower dimensional element corresponding to
      // the secondary face. This must be done last after the dof indices have been prepared for the
      // secondary (element) and primary (neighbor)
      subproblem.reinitLowerDElem(secondary_face_elem, /*tid=*/tid, &xi1_pts);

      // All this does currently is sets the neighbor/primary lower dimensional elem in Assembly and
      // computes its volume for potential use in the MortarConstraints. Solution continuity
      // stabilization for example relies on being able to access the volume
      subproblem.reinitNeighborLowerDElem(msinfo.primary_elem, tid);

      // reinit higher-dimensional secondary face/boundary materials. Do this after we reinit
      // lower-d variables in case we want to pull the lower-d variable values into the secondary
      // face/boundary materials. Be careful not to execute stateful materials since conceptually
      // they don't make sense with mortar (they're not interpolary)
      fe_problem.reinitMaterialsFace(secondary_ip->subdomain_id(),
                                     /*tid=*/tid,
                                     /*swap_stateful=*/false,
                                     &secondary_ip_mats);
      fe_problem.reinitMaterialsBoundary(
          secondary_boundary_id, /*tid=*/tid, /*swap_stateful=*/false, &secondary_boundary_mats);

[Jensoulou]

I'm sorry for not getting back to you sooner. I'm still trying to figure it all out.

My goal is to adapt the mechanical PorousFlow kernels, such as PorousFlowMassVolumetricExpansion and PorousFlowHeatVolumetricExpansion, to work on LowerD elements. If I simply apply them as they are now on a LowerD subdomain, the computation is wrong: these kernels rely on the strain rate and its derivative, which in turn require the full-dimensional strain. In LowerD elements, one direction is missing, so the strain rate is computed incorrectly compared to the HigherD elements at the same location. Physically, these values should be identical across the two representations.

On top of this adaptation, I also want to introduce an AuxVariable for the aperture of the fault modeled with LowerD elements. This AuxVariable would feed into the adapted mechanical PorousFlow kernels, allowing me to modulate the fault’s influence based on its aperture, which in the future would evolve with stress/strain or other variables.

Thank you for the pointer to this snippet from MortarUtils.h. From what I understand, the Lower-D elements are reinitialized after the Higher-D elements and their faces, which ensures that the data related to the Higher-D side is already up to date. I am, however, not fully certain about the following, and would like to check with you if my understanding is correct.

We stated earlier that the Lower D and Higher D elements are not neighbors in the mesh sense, and cannot be treated as such. But, in the mortar code, since both are set up in the same routine (fe_problem.reinitNeighborFaceRef(...) then subproblem.reinitLowerDElem(...)), both contexts have been prepared when we are looking from the Lower D's point of view. This means that the Higher-D data has already been stored in the Assembly. This data can be accessed via
the neighbor functions, such as getNeighborMaterialProperty<Real>("strain_rate"), because loopOverMortarSegments routine establishes a virtual link that treats the Higher-D element face as the “neighbor” of the Lower-D element. So, although they are not neighbors from the mesh topology perspective, the mortar code uses the neighbor functions to allow the Lower-D element to access data from the Higher-D element.

Does this interpretation align with what's happening in the mortar code? If so, I believe the mortar code could perform the task I need.

[lindsayad]

You're definitely doing a good job reading the code! But that interpretation isn't quite right so good to check. We don't establish a neighbor connection because as you point out that would be a bit silly from a topological standpoint. It's just that we are using the same physical space quadrature points for reinitializing both the higher-dimensional element and lower-dimensional element so data can then be used from both. Perhaps the ADMortarConstraint constructor could be helpful to see

ADMortarConstraint::ADMortarConstraint(const InputParameters & parameters)
  : MortarConstraintBase(parameters),
    _lambda_dummy(),
    _lambda(_var ? _var->adSlnLower() : _lambda_dummy),
    _u_secondary(_secondary_var.adSln()),
    _u_primary(_primary_var.adSlnNeighbor()),
    _grad_u_secondary(_secondary_var.adGradSln()),
    _grad_u_primary(_primary_var.adGradSlnNeighbor())
{
  _subproblem.haveADObjects(true);
}

_lambda is a variable that lives only on the lower-dimensional elements. We retrieve its solution at quadrature points using the adSlnLower() variable API. (It would be simply slnLower() if we were hand-forming Jacobians). And then the primal solution on the higher-dimensional secondary face (which is the face coincident with lower-d block) is retrieved using the adSln() API. (The mortared side (the other higher-dimensional face we are pairing) is the topological neighbor so we retrieve with adSlnNeighbor() there).

Both adSlnLower() and adSln() will be of the same length--the number of mortar quadrature points--after the sequence of reinit calls I pasted previously. So they can then be used in the same residual computing code paths. You can see

moose/framework/src/constraints/EqualValueConstraint.C

Lines 43 to 97 in 17a5834

	ADReal
	EqualValueConstraint::computeQpResidual(Moose::MortarType mortar_type)
	{
	switch (mortar_type)
	{
	case Moose::MortarType::Secondary:
	{
	// The sign choice here makes it so that for the true solution: lambda = normals_secondary *
	// diff_secondary * grad_u_secondary
	auto residual = -_lambda[_qp] * _test_secondary[_i][_qp];
	if (_stabilize)
	residual += _delta * _lower_secondary_volume *
	(_diff_secondary[_qp] * _grad_test_secondary[_i][_qp] * _normals[_qp]) *
	(_lambda[_qp] - _diff_secondary[_qp] * _grad_u_secondary[_qp] * _normals[_qp]);
	return residual;
	}
	case Moose::MortarType::Primary:
	{
	// The sign choice here makes it so that for the true solution: lambda = -normals_primary *
	// diff_primary * grad_u_primary
	auto residual = _lambda[_qp] * _test_primary[_i][_qp];
	if (_stabilize)
	residual -= _delta * _lower_primary_volume *
	(_diff_primary[_qp] * _grad_test_primary[_i][_qp] * _normals[_qp]) *
	(_diff_primary[_qp] * _grad_u_primary[_qp] * _normals[_qp] - _lambda[_qp]);
	return residual;
	}
	case Moose::MortarType::Lower:
	{
	auto residual = (_u_primary[_qp] - _u_secondary[_qp]) * _test[_i][_qp];
	if (_stabilize)
	{
	// secondary
	residual -= _delta * _lower_secondary_volume * _test[_i][_qp] *
	(_lambda[_qp] - _diff_secondary[_qp] * _grad_u_secondary[_qp] * _normals[_qp]);
	// primary
	residual -= _delta * _lower_primary_volume * _test[_i][_qp] *
	(_lambda[_qp] - _diff_primary[_qp] * _grad_u_primary[_qp] * _normals[_qp]);
	}
	return residual;
	}
	default:
	return 0;
	}
	}

for an example