OpenCMISS / iron

Source code repository for OpenCMISS-Iron
9 stars 62 forks source link

Neumann_point BC type #128

Open beauof opened 6 years ago

beauof commented 6 years ago

The Neumann_point BC type gives unexpected results.

For a simple uniaxial extension test, the results look as if one does not apply consistent nodal forces, however, that's what this BC type is supposed to compute.

Needs debugging or proper documentation on how to use this BC type!

chrispbradley commented 6 years ago

Hi Andreas, do you have an example of zip file that you could attach that shows the problem. What RHS values are you expecting?

beauof commented 6 years ago

The example that I am currently working on is this one: https://github.com/beauof/iron-tests/blob/master/examples/example-0201-u/src/iron/example.F90 and it is documented (sparsely at the moment) here: https://github.com/beauof/iron-tests/blob/master/tests.pdf

It's for a 3D uniaxial test:

It uses the following variables to set up the simulation:

! defaults for input arguments WIDTH = 1.0_CMISSRP HEIGHT = 0.2_CMISSRP LENGTH = 0.2_CMISSRP NumberGlobalXElements = 2 NumberGlobalYElements = 2 NumberGlobalZElements = 2 SolverIsDirect = 1 ! direct solver by default JACOBIAN_FD = 1 ! finite-difference Jacobian by default MooneyRivlin1 = 0.5_CMISSRP*35.7_CMISSRP ! see note above for Neo-Hookean solid MooneyRivlin2 = 0.0_CMISSRP ! If MooneyRivlin2 == 0 --> Neo-Hookean solid useGeneratedMesh = 1 ! generated mesh by default BCDISP_MAX = 0.2_CMISSRP bcType = 0 ! 0 - Dirichlet BC by default; else: 1 - Neumann_integrated, 2 - Neumann_point

A BCDISP_MAX=0.2 corresponds to a stretch of 1.2.

A bcType=0 will run the displacement-driven problem. A bcType=1 will run the problem with Neumann_integrated BC and bcType=2 is meant to use the Neumann_point BC.

There are a few routines that read user-defined meshes if requested, that you can comment out. There is also a routine that computes the nodal weights for the Neumann_integrated BC necessary to run bcType=1. The routines are here at the moment: https://github.com/beauof/meshReader/blob/master/src/meshReader.F90 But you won't need them if you run the problem with generated meshes and bcType=1 or bcType=2.

Thanks for having a look!

lorenzo-mechbau commented 5 years ago

For this type of traction-BC the trick is to use Neumann_point on the traction surface (with input force values given in the current configuration) AND apply 0. Neumann_integrated_only conditions everywhere else. This way, integration happens only on the required Neumann surface, otherwise it occurs on ALL faces adjacent to every Neumann node.

The way Neumann_point has been coded is meant to deal with forces applied in one point, but affecting the whole surface. In case of distributed forces you can either use this workaround or re-code the Neumann integration (I have a branch for that).