Project IP to Nodal Field fails with Composite Tets

[lxmota]

@ikalash

Project IP to Nodal Field for computing an L2 projection of the stress to the nodes fails for composite tetrahedral elements.

This is a long-standing issue that we need to fix, as users keep running into it.

[ikalash]

Just to clarify: my concern / frustration was that an error is not thrown when trying to use "Project IP to Nodal Field" with composite tets. I got a generic error from Intrepid2 about the cell topology not being supported and thought for awhile that there was something messed up with my tet10 mesh. If there is no plan to extend "Project IP to Nodal Field" to composite tets, that's fine, but the code should throw an error that this response is not supported when trying to use it with composite tets.

[mperego]

If you need the projection to work with composite tets and you think it is an issue with Intrepid2 let me know and I can look into it

[jtostie]

Fixing the underlying issue in the composite tet projection involves using a higher order numerical integration scheme to compute the projection matrix. I cannot commit to work on that. In another code we exploit the assumed linear character of the integration point field and project using the linear tet4 basis.

Once upon a time there was another response option, IP to Nodal Field, (without the Project) that worked consistently with the composite tet10's in Albany. That may serve as an alternative.

Given that, it seems like we should just disallow the Project* response when using tet10s.

[lxmota]

@jtostie That's OK, I just assigned to you and me because I remember having discussions about this.

The IP to Nodal Field option works, so I guess that for now I'll go ahead and disable Project IP to Nodal Field for composite tets.

[ikalash]

A follow up comment / question for @mperego . The error that is thrown when we try to use the Project IP to Nodal Field response with composite tets is the following.

Intrepid2] Error in file /home/ikalash/Trilinos/build-dtk-release/install/include/Intrepid2_CellTools.hpp, line 200
Test that evaluated to true: true
>>> ERROR (Intrepid2::CellTools::createHGradBasis): Cell topology not supported.

p=0: *** Caught standard std::exception of type 'std::invalid_argument' :

ERROR (Intrepid2::CellTools::createHGradBasis): Cell topology not supported.

Looking at Intrepid2_CellTools.hpp, I am wondering, why is line 180 commented out:

//case shards::Tetrahedron<11>::key: r_val = Teuchos::rcp(new Basis_HGRAD_TET_COMP12_FEM<ExecSpaceType,outputValueType,pointValueType>()); break;

I think the element with the tet-11 key is the composite tet one and this line being commented out might be the problem.

[ambrad]

I stumbled across this issue and am copy-(cleaning-)pasting some notes I made on
this topic way back in June 2015, in case they are of any use:

1. For a linear test function, tet10 recovers the test function exactly.

2. If the tet10s are straight-edged with edge nodes at the middle of the edges,
   then tet10 can and does recover a quadratic test function exactly.

3. But if the edge node is not at the center, then it can't. The Jacobian is a
   rational function and so exactness can't be expected.

4. Increasing the quadrature order should fix this, but something made me
   conclude that quadrature for composite tet10 is not supported beyond degree 4,
   or something roughly like that.

5. Last time I worked on this, I don't recall there being an Intrepid issue. So
   I'm guessing there's just a currently unconverted part of Intrepid2, or
   something like that.

[ikalash]

@mperego : Per our discussion at the Albany meeting this week, we definitely do need support for the LCM composite tet element in Intrepid2 (tagging @lxmota who can comment more if necessary). I have looked in the Intrepid2 source code and it looks like most of the composite tet stuff is commented out. Can you or @kyungjoo-kim please port it there from Intrepid? Since the Albany::ProblemUtils class constructs its own Intrepid2 basis with the relevant basis functions, we can run LCM problems with composite tets; but then when we try to do some analysis that requires other Intrepid2 functionality specific to the composite tets (e.g., the Project IP to Nodal Field), we encounter issues.

[mperego]

It seems that the code for composite tets has been ported to Intrepid2 and tested (at least partially), except for the quadrature rule CubatureCompositeTet, that it is not used anyway (even in Intrepid).
Can somebody that knows what to expect try to uncomment these two lines:
https://github.com/trilinos/Trilinos/blob/master/packages/intrepid2/src/Cell/Intrepid2_CellTools.hpp#L80
and
https://github.com/trilinos/Trilinos/blob/master/packages/intrepid2/src/Cell/Intrepid2_CellTools.hpp#L180
and see what happens?

I'm on travel next week but I can dig into this after that.

[ikalash]

@mperego: when I uncomment lines 80, 139 and 180 in Intrepid2_CellTools.hpp and run an Albany example with the composite tet element + the Project IP to Nodal Field Response, the problem runs; in contrast to before where we got a cryptic seg fault. So I think this uncommented code would be used with composite tet elements + Project IP to Nodal Field response, but there are probably no examples in Albany that test this (because this response is not accurate with the given cubatures, as discussed above). I would advise you to check in this change to Intrepid2 to make Intrepid2 consistent with Intrepid, to avoid cryptic seg faults in Albany, and (maybe, if it is possible) to have a place to have a place to add new cubatures for the composite tet in the future.

[ambrad]

@ikalash, if it's of use (and still works), I believe every output of ProjectIPToNodalField includes a linear test function with "linear" in the Exodus field name. This can be used to sanity check results, since a linear function should I think in all reasonable circumstances be recovered exactly.

[ikalash]

Thanks, @ambrad . I actually found this while looking through the code within the PROJ_INTERP_TEST ifdef guards.

[mperego]

@ikalash It is now possible to select the composite Tet basis using Intrepid2::CellTools::createHGradBasis method (see trilinos/Trilinos@f06aca8)

[ikalash]

OK, thanks for propagating that code to Intrepid2, @mperego !