JxWNeighbor() method make interfacekernel crash

[arovinelli]

Bug Description

Assembly's method JxWNeighbor cannot be used to get neighboring JxW. I assume something is not initialized or reinited properley @lindsayad

Steps to Reproduce

In the InterfaceKernel.C constructor add _JxW_neighbor(_assembly.JxWNeighbor()) In the InterfaceKernel.C::computeElemNeighResidual add

  const MooseArray<Real> &loc_JxW =
      (type == Moose::Element) ? _JxW : _JxW_neighbor;

  for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     std::cout << "_qp" << _qp << "_JxW_neighbor " << _JxW_neighbor[_qp];

and it will crash with a stacktrace similar to

Assertion `i < _size' failed
Access out of bounds in MooseArray (i: 0 size: 0)
at /home/arovinelli/projects/NEW_stuff/moose/framework/build/header_symlinks/MooseArray.h, line 277
Detaching after fork from child process 36659.
Stack frames: 26
0: libMesh::print_trace(std::ostream&)
1: MooseArray<double>::operator[](unsigned int) const
2: InterfaceKernelTemplDisconnectedMesh<double>::computeElemNeighResidual(Moose::DGResidualType)
3: InterfaceKernelTemplDisconnectedMesh<double>::computeResidual()
4: ComputeResidualThread::onInterface(libMesh::Elem const*, unsigned int, short)
5: ThreadedElementLoopBase<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> >::operator()(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&, bool)
6: void libMesh::Threads::parallel_reduce<libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*>, ComputeResidualThread>(libMesh::StoredRange<libMesh::MeshBase::const_element_iterator, libMesh::Elem const*> const&, ComputeResidualThread&)
7: NonlinearSystemBase::computeResidualInternal(std::__debug::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&)
8: NonlinearSystemBase::computeResidualTags(std::__debug::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&)
9: FEProblemBase::computeResidualTags(std::__debug::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&)
10: FEProblemBase::computeResidualInternal(libMesh::NumericVector<double> const&, libMesh::NumericVector<double>&, std::__debug::set<unsigned int, std::less<unsigned int>, std::allocator<unsigned int> > const&)

Impact

It does prevent me to finalize a DG based large deformation CZM model

[lindsayad]

What's the idea here? That interfacial areas are different between slave and master elements? With this CZM stuff in general the physical location of the integration points for the element are not actually guaranteed to lie on the neighbor are they? If that's the case I'm not actually sure how any of this CZM stuff works because of this code:

void
Assembly::reinitElemAndNeighbor(const Elem * elem,
                                unsigned int side,
                                const Elem * neighbor,
                                unsigned int neighbor_side)
{
  _current_neighbor_side = neighbor_side;

  reinit(elem, side);

  unsigned int neighbor_dim = neighbor->dim();

  std::vector<Point> reference_points;
  FEInterface::inverse_map(
      neighbor_dim, FEType(), neighbor, _current_q_points_face.stdVector(), reference_points);

  reinitFEFaceNeighbor(neighbor, reference_points);
  reinitNeighbor(neighbor, reference_points);
}

where we try to determine what the reference points are on the neighbor based on the physical location of the integration points coming from the elem.

[lindsayad]

I'm anticipating error messages out of here:

      // Make sure the point \p p on the reference element actually
      // is

      if (!FEAbstract::on_reference_element(p, elem->type(), 2*tolerance))
        {
          libmesh_here();
          libMesh::err << "WARNING:  inverse_map of physical point "
                       << physical_point
                       << " is not on element." << '\n';
          elem->print_info(libMesh::err);
        }

when running in debug mode. Do you ever see this in debug mode? Or are your slave and master surfaces always overlapping?

[arovinelli]

What's the idea here? That interfacial areas are different between slave and master elements?

yes they are different, up to now we didn't consider such difference but area changes on both side as deformation evolves.

With this CZM stuff in general the physical location of the integration points for the element are not actually guaranteed to lie on the neighbor are they?

They are not guaranteed to lie on top of each other. The are coincident only when disaplcement are 0

where we try to determine what the reference points are on the neighbor based on the physical location of the integration points coming from the elem.

I remeber @roystgnr modified some code in libmesh for the same reason. It was making fail the breakmeshblock generator in debug mode libmesh/pull/2363

I'm anticipating error messages out of here:

not yet :)

[srinath-chakravarthy]

I am getting these messages in the debug mode ?

Andrea, Could you share the reference you are using for the large deformation czm.

Cheers Srinath

On Thu, Mar 12, 2020, 4:13 PM Alex Lindsay notifications@github.com wrote:

I'm anticipating error messages out of here:

  // Make sure the point \p p on the reference element actually
  // is

  if (!FEAbstract::on_reference_element(p, elem->type(), 2*tolerance))
    {
      libmesh_here();
      libMesh::err << "WARNING:  inverse_map of physical point "
                   << physical_point
                   << " is not on element." << '\n';
      elem->print_info(libMesh::err);
    }

when running in debug mode. Do you ever see this in debug mode? Or are your slave and master surfaces always overlapping?

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/idaholab/moose/issues/14870#issuecomment-598394710, or unsubscribe https://github.com/notifications/unsubscribe-auth/ACYC6LJDXDP3NCNSTDAPX6LRHE66JANCNFSM4LFI2BFQ .

[arovinelli]

I am getting these messages in the debug mode ? Andrea, Could you share the reference you are using for the large deformation czm. Cheers Srinath

Dear Srinath @srinath-chakravarthy, I'm developing the CZM large deformation model in these days and I'm still in the process of working out some issues with the MOOSE team. Once I have a pre-release working version I'll let you know.

[lindsayad]

So the fact that _JxW_neighbor makes interface kernels crash is maybe a bug (I'll have to see whether there's a way to error up front if a developer tries to use JxwNeighbor in a DGKernel or InterfaceKernel), but you should never be using _JxW_neighbor in interface kernels to begin with. That member exists for use in things like NodeFaceConstraint (see its original creation here: https://github.com/idaholab/moose/issues/3317). By trying to use sides which are not coincident, you are violating the concept of the interface kernel. We are back to talking about using objects based on geometric searches

[arovinelli]

So the fact that _JxW_neighbor makes interface kernels crash is maybe a bug (I'll have to see whether there's a way to error up front if a developer tries to use JxwNeighbor in a DGKernel or InterfaceKernel), but you should never be using _JxW_neighbor in interface kernels to begin with. That member exists for use in things like NodeFaceConstraint (see its original creation here: #3317). By trying to use sides which are not coincident, you are violating the concept of the interface kernel. We are back to talking about using objects based on geometric searches

@lindsayad I understand your point but I don't think #14890 is the right way to go. This problem could be solved with a geometric search (as you mentioned). DG with disconnected elements it's use frequently in literature for large displacement disconnected mesh problems. Do you think it is reasonable for me to code a CohesiveGeomtricSearch that identifies neighbor at the beginning and then simply returns that at every time step. It shouldn't be some much overhead (only computed once, then it's stored somewhere, like UO) and then probed anytime you want to do you your costly inversion?
I see in code MOOSE already has a ElementPairInfo and mortar constraint that with some work could solve all this issues, allowing CZM to work with InterfaceKernel and displaced mesh. In this case, mortar would be more a fixed neighbor.

[lindsayad]

For the interfacial physics I’ve worked with involving mesh displacement (thermal contact and mechanical contact), you always want to work with your geometrically closest neighbor. Intuitively this makes sense to me.

For cohesive zone modeling, you really use the original neighbor regardless of the amount of mesh deformation? If that’s the case then we can work something out that’s fairly close to DGKernel/InterfaceKernel

On Mar 14, 2020, at 4:58 PM, Andrea Rovinelli notifications@github.com wrote:

 So the fact that _JxW_neighbor makes interface kernels crash is maybe a bug (I'll have to see whether there's a way to error up front if a developer tries to use JxwNeighbor in a DGKernel or InterfaceKernel), but you should never be using _JxW_neighbor in interface kernels to begin with. That member exists for use in things like NodeFaceConstraint (see its original creation here: #3317). By trying to use sides which are not coincident, you are violating the concept of the interface kernel. We are back to talking about using objects based on geometric searches

@lindsayad I understand your point but I don't think #14890 is the right way to go. This problem could be solved with a geometric search (as you mentioned). DG with disconnected elements it's use frequently in literature for large displacement disconnected mesh problems. Do you think it is reasonable for me to code a CohesiveGeomtricSearch that identifies neighbor at the beginning and then simply returns that at every time step. It shouldn't be some much overhead (only computed once, then it's stored somewhere, like UO) and then probed anytime you want to do you your costly inversion?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[arovinelli]

For cohesive zone modeling, you really use the original neighbor regardless of the amount of mesh deformation? If that’s the case then we can work something out that’s fairly close to DGKernel/InterfaceKernel

@lindsayad

It really depends on what you are trying to model, and basically there are two cases:

1. elastic interface medium, upon unloading everything goes back to the initial configuration, then the physical neighbors should never change even in presence of large deformations

2. visco elastic medium, in which sliding is prominent, than physical neighbors should change. However, for cohesive modeling we need to avoid what is the master slave concepts present in mortar or contact, for which we do need a refine mesh on one sided and a coarse mesh on the other.

The best would be to have a single kernel that would allow to switch between the two behaviors. My best guess now would be to use a facefaceconstraint, in which we update physical neighbors only for case 2. For now I would be more than happy to have case 1 working properly. I know you are busy, but it would be very beneficial if we could setup some time to talk and setup a pathway that is optimal in the short and long term for INL, ANL and MOOSE. It shouldn't take long.

[lindsayad]

If cohesive zone is being setup with an action, which I believe it currently is, it doesn’t really matter whether 1 and 2 are done by the same object. You could toggle between the two with an action parameter and then the action could add a different object if say an interfacekernel is the more natural treatment for 1 while a mortar constraint or something akin to it would be a more natural treatment for 2. Given that, let’s go ahead with 1. I think I can have a solution done up for that in a couple hours (when I get to it, probably sometime this coming week).

On Mar 15, 2020, at 10:49 AM, Andrea Rovinelli notifications@github.com wrote:

 For cohesive zone modeling, you really use the original neighbor regardless of the amount of mesh deformation? If that’s the case then we can work something out that’s fairly close to DGKernel/InterfaceKernel

@lindsayad

It really depends on what you are trying to model, and basically there are two cases:

elastic interface medium, upon unloading everything goes back to the initial configuration, then the physical neighbors should never change even in presence of large deformations

visco elastic medium, in which sliding is prominent, than physical neighbors should change. However, for cohesive modeling we need to avoid what is the master slave concepts present in mortar or contact, for which we do need a refine mesh on one sided and a coarse mesh on the other.

The best would be to have a single kernel that would allow to switch between the two behaviors. My best guess now would be to use a facefaceconstraint, in which we update physical neighbors only for case 2. For now I would be more than happy to have case 1 working properly. I know you are busy, but it would be very beneficial if we could setup some time to talk and setup a pathway that is optimal in the short and long term for INL, ANL and MOOSE. It shouldn't take long.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[arovinelli]

If cohesive zone is being setup with an action, which I believe it currently is, it doesn’t really matter whether 1 and 2 are done by the same object. You could toggle between the two with an action parameter and then the action could add a different object if say an interfacekernel is the more natural treatment for 1 while a mortar constraint or something akin to it would be a more natural treatment for 2. Given that, let’s go ahead with 1.

Yes we use an action so we can definitely switch between the two giving some flag. Yes you are right, for sliding we definitely need to think something else. Let's put 2 on hold for now.

I think I can have a solution done up for that in a couple hours (when I get to it, probably sometime this coming week).

This would be great. I already have a material doing all the kinematics and as soon as you can fix this I can make a PR to enable large deformation czm without sliding. This will be particularly useful for people crack propagation and stuff like that, where sliding is irrelevant

[jiangwen84]

@arovinelli I had a discussion with Ben and Derek about a year ago regarding to CZM large deformation and/or large sliding. We think FaceFaceConstraint would be better suited for that. There are more than the two cases that you listed above. With FaceFaceConstraint, you should be able to handle all the cases within a single framework.

[lindsayad]

Implementing the first scenario @arovinelli outlined would be quite a bit slower using anything that would require geometric searches. While something like FaceFaceConstraint would be a lot more general and probably more physically accurate, if the modeler is happy with a tied neighbor approach, I think the computational science team can honor their request