Higher order integration for Neumann BC

[andrea-franceschini]

What is the requested feature? Currently, we impose traction boundary condition as F*A/n, with n number of nodes for the specific face. We should use a proper integration. E.g.: https://github.com/GEOSX/GEOSX/blob/ee9de4c0693462b2471567882aeb62f16d452c1c/src/coreComponents/physicsSolvers/solidMechanics/SolidMechanicsLagrangianFEM.cpp#L736

Is your request related to a specific problem? For quadrilateral faces with generic angles, this is a rough approximation.

Describe the solution you'd like To use a proper integration, from the new finite element library.

[castelletto1]

We should address this as part of #1094