opencmiss / iron Goto Github PK

This issue tracks cleanup and incorporation of my remaining PhD code into iron:develop. This includes 3D-0D coupling and fixes for boundary face integration, as well as any relevant bugs that turn up.

The branch (dladd/iron:dladd_phd) that this work is in was pulled/merged then reverted so it may eventually need to be manually merged (see discussion in issue #74). But I'll try to update this parent issue to track my progress before attempting to re-merge.

TODO:

• Boundary element line/face integration
  • Generalise boundary line/face integration to 2D/3D (currently 3D faces only)
  • Generalise for simplex element types
• Remove legacy/nonfunctioning 3D-1D code
• Re-work/generalise hard-coded fix calling CellML solver at each residual evaluation (see changes to problem_routines.)
• Fix issue #76 : a bug in Navier-Stokes code expecting a specific ordering of the velocity and pressure mesh components
• Change implementation of pressure boundary conditions to a proper Neumann type rather than going through the PRESSURE_VALUES parameter set.

Iron code needs to be converted to the new CaMeL case code style.

Increase range of models that can be read and written by Iron in FieldML 0.5 format.

The equations and solver equations should "register" what parameter sets are required as part of the mapping setup. The mapping would then collate the parameter sets for the field variables. This way the parameter sets can be ensured to be created when the mapping is finished.

It should be possible to destroy a mesh once it has been decomposed. This will save a fair bit of memory. There are, however, a few references to mesh attributes. Before the mesh can be destroyed these references need to be changed to references to decomposition attributes. The decomposition may also have to be moved from a mesh to a region as the meshes currently hold the list of decompositions.

The essential message is PyThreadState_Get: no current thread
As the other tests regarding Python pass, i think that has to do with the virtual environment there. Needs checking.

59/59 Test #59: Python_Bindings_Virtualenv ...............................***Failed 23.09 sec
New python executable in /Users/danielthewicked/opencmiss/build/x86_64_darwin/clang-6.0-F4.5/mpich_release/iron/debug/bindings/python/tests/test_venv/bin/python2.7
Also creating executable in /Users/danielthewicked/opencmiss/build/x86_64_darwin/clang-6.0-F4.5/mpich_release/iron/debug/bindings/python/tests/test_venv/bin/python
Installing setuptools, pip, wheel...done.
Processing /Users/danielthewicked/opencmiss/build/x86_64_darwin/clang-6.0-F4.5/mpich_release/iron/debug/bindings/Debug
Installing collected packages: OpenCMISS-Iron
Running setup.py install for OpenCMISS-Iron: started
Running setup.py install for OpenCMISS-Iron: finished with status 'done'
Successfully installed OpenCMISS-Iron-0.5.0
Fatal Python error: PyThreadState_Get: no current thread
CMake Error at /Users/danielthewicked/opencmiss/src/iron/bindings/python/tests/TestVirtualenv.cmake:26 (message):
Importing opencmiss iron package in virtual environment failed

Currently when fields are exported to exfiles no validation is performed on iron field variable names to check if they meet requirements of zinc. This is an issue for me with the "del U/del n" field ending up with a field name containing the region path separator.

A simple solution is to simply replace '/' with '_' when exporting. Which at least results in a field name that can be used in zinc (actually, probably this is a Cmgui only issue?).

cmiss.h is defining SA_NODEFER but on OS X 10.11 the system has already done so

We have now moved from a triplet of integers (class, type and subtype) to an arbitrary length integer array for specifying equations sets and problems. We are now thus in a position to start to add more integers in order to clean up some of the specifications which have overloaded subtype with a vast number of different options. A good example of the problem is mechanics. The problem is how does the user know how big an integer array is required? Do they have to specify all options in the array? What are the defaults? What if another option/array slot is required in the future? In order to allow for future expansion I propose that we have a specification object that is setup via a create start, create finish mechanism. The create start would allocate the specification array (inside the specification object) and set up all defaults. Different options for the specification would then be set by giving an integer parameter that indicates the option slot and another integer parameter for the option value. Create finish would then finish the object i.e.,

EquationsSetSpecification = iron.EquationsSetSpecification()
EquationsSetSpecification.CreateStart(EquationsSetSpecificationUserNumber)
EquationsSetSpecification.OptionSet(iron.EquationsSetOptionType.CLASS,iron.EquationsSetOptionValue.ELASTICITY_CLASS)
EquationsSetSpecification.OptionSet(iron.EquationsSetOptionType.TYPE,iron.EquationsSetOptionValue.ELASTICITY_TYPE)
EqautionsSetSpecification.OptionSet(iron.EquationsSetOptionType.CONSTITUATIVE_LAW, iron.EquationsSetOptionValue.MOONEY_RIVILIN_CONSTITUATIVE_LAW)
EquationsSetSpecification.OptionSet(iron.EquationsSetOptionType.FINITE_ELASTICITY_GROWTH_LAW,iron.EquationsSetOptionValue.NO_GROWTH_LAW)
EquationsSetSpecification.CreateFinish()

etc.

The slight issue with this approach is that the option types and constants could not be "enums" in that multiple sets would have to be used for the same argument slot in the call. A giant list of identifies would thus be required. It may be that doxygen could be used to group the constants into more sensible groups?

The same approach would be used for problem specification.

It might be better to move to a relative specification of field variables in equations sets. At the moment field variables are referenced by their variable type in an equation. In the relative mapping an equation knows it has n variables and just references them as variable 1, variable 2 etc. The mapping then takes care of what actual field variable is used. This will allow more flexibility for coupled and multi-physics problems.

Iterative solvers diverge for the generalised Laplace problem in 2D.

The error is caused by uninitialised variables in the finite element calculate routines, a bug fix is available (see pull-request #82 )

I ran into a bug (segfault) when using DECOMPOSITION_TOPOLOGY_FACES_CALCULATE with simplex element types. The relevant section of code is in mesh_routines around L3200:

DECOMPOSITION_FACE%XI_DIRECTION=BASIS%LOCAL_FACE_XI_DIRECTION(basis_local_face_idx)
DOMAIN_FACE%BASIS=>BASIS%FACE_BASES(ABS(DECOMPOSITION_FACE%XI_DIRECTION))%PTR
ALLOCATE(DOMAIN_FACE%NODES_IN_FACE(BASIS%NUMBER_OF_NODES_IN_LOCAL_FACE(basis_local_face_idx)), &
& STAT=ERR)
IF(ERR/=0) CALL FlagError("Could not allocate face nodes in face",ERR,ERROR,*999)
ALLOCATE(DOMAIN_FACE%DERIVATIVES_IN_FACE(2,DOMAIN_FACE%BASIS%MAXIMUM_NUMBER_OF_DERIVATIVES, &
& BASIS%NUMBER_OF_NODES_IN_LOCAL_FACE(basis_local_face_idx)),STAT=ERR)

For a tetrahedron, we can't allocate DOMAIN_FACE%BASIS since BASIS%FACE_BASES range from 1:3 but DECOMPOSITION_FACE%XI_DIRECTION range from 1:4

I mentioned this at the AKL meeting and it was thought a likely cause of this may be confusion of XI vs. XIC directions for simplex types.

cases for "irix" and "vms" are not yet covered (machine_constants_XXX.f90) within the cmake build system

First part in the effort of reducing the physical memory requirement in Iron. Want to remove the use of the GLOBAL_TO_LOCAL_MAP structure in all domain mappings (eg. element, node and DOF mappings). Idea is to use the LOCAL_TO_GLOBAL_MAP array instead for global number look-ups.

The finite elasticity code is a bit of a dog's breakfast at the moment. The code needs to be tidied and reorganised. Some steps include

• Sort out a hierarchy of equations set specifications to capture all the various optons.
• Allow a more modular nature to the the elasticity stress calculation e.g., allow an arbitrary number of sub-equations types which can each calculate part of the Cauchy stress. The final Cauchy stress is then a sum of all the sub-parts. Some sub-parts could include elastic, visco-elastic, active, passive etc. etc.
• Also need to introduce a temporal aspect in order to accommodate viscoelasticity, plasticity, dynamics etc.

For a transient NS problem it is common to have quadratic for velocity and linear for pressure. Currently NS works if the velocity (quadratic) mesh component is number 1 and the pressure (linear) mesh component is number 2. The NS code does not work, however, if you have the linear pressure component as mesh component 1 and the quadratic velocity component is mesh component 2. There is obviously some assumptions going on in the code which needs to be fixed.

This issue is intended for tracking progress with respect to reproducing the tibialis anterior model published in Heidlauf et al. 2013.

We (@beauof, @harnoorsaini, @klotz-t) will start with commit ea441e5 and add details to this issue once we identified problems.

• fix compilation error
• track down bugs

edited by @beauof for tracking progress:

TA model does not work with:

• commit ea441e5
• Prasad's finite-elasticity-fixes branch
• ...

In order to be able to really develop C examples against the iron C bindings, the correct dllimport declarations have to be used. CMake offers GenerateExportHeader etc for that.
it just has not been done yet.

At the moment a region can only hold one set of data points. This needs to be extended to allow for multiple sets of data points. This will mean that data points will need to have a user number within a region.

• At the moment there is only one BC flag array. This was meant to be just for the DOF being FIXED or FREE but it has now become overloaded with various different BC 'types'. We need another array for the BC 'type' which then fixes the FIXED/FREE for the particular DOF.

It might be better to have a way of creating field hierarchies. Field variables are made up of a number of components. For some problems these components are mixed e.g., for mechanics problems we have displacement/velocity and pressure bundled together. It would be better to have a hierarchy where we explicitly group the displacement/velocity "field" with the pressure "field". This would allow the displacement/velocity field to keep its geometrical vector property. There should be an arbitrary nature to the bundling i.e., the fields are bundled with components (field variables) which are then bundled with other fields bundles of components (field) which are then bundled (fields) etc. This would then move the world/object field to that more matching the mathematical definition of scalar, vector, tensor fields etc.

Need a convenient routine to (re)normalise derivatives so that they have unit length. An issue that may have to be addressed first is that each field component can use a different mesh component and so there is no reason to assume that a node will exist in all geometric field components or that it is even the same local node number. It may be that we need to store in a decomposition the list of "mesh" nodes and what mesh components that node is in.

Need to add the public APIs for setting/getting the newton solver solution init type.

Issue for the discussion of what Fortran 2003 features can be used in OpenCMISS. In general, the use of language features depends on performance implications and the support of that feature by the various Fortran compilers. A useful summary of Fortran language features that have been implemented by the various compilers can be found here:

http://fortranwiki.org/fortran/show/Fortran+2003+status

To fully attain model interchange between Iron and Zinc we need to head towards a common or at least compatible data model/s for element parameter mapping and scale factor indexing. At the same time we should implement general linear mapping which enables elegant methods to deal with apexes, hanging nodes, managing continuity and more.

Following is a link to a proposed data model that we can head towards; it describes the aims, current status of element parameter mapping, requirements, examples of mappings including general linear maps, and the proposal itself:

https://docs.google.com/document/d/1giu0nNSIm2prg4yQ9G3FzkDfzQ6_jFE82vOnt7BneI4/edit?usp=sharing

We'd appreciate your feedback. Questions/comments and suggestions for improving the text can be put in the document. Your overall indication of support, or major concerns and alternative ideas can be added to this issue.

When an iron error is encountered with python bindings the response is to raise an exception and output the error string obtained with ExtractErrorMessage. This is slightly different to other bindings which call cmfe_HandleError. One consequence is that the stack trace is not output and thus there is no information on what routine that error occurred in. The error behaviour should be a bit more consistent and that stack trace should be output for python.

I was wondering what the guideline is regarding reading third party meshes, for example from Cubit/Abaqus/etc. file format.

Recently, the developers in Stuttgart have written custom mesh routines for reading CHeart meshes (one of the formats that we frequently use) and it would be great to include that functionality in OpenCMISS-iron and in opencmiss_iron.f90 specifically, rather than copy'n'pasting those routines to every single example file.

We would like to avoid
(i) converting and thus duplicating meshes before reading into iron and
(ii) copy'n'pasting routines to new example files.

Any comments/suggestions?

Check the makefile-based compiler flags file to check and import the correct flags for AIX systems

DataPointsProjectionEvaluate segfaults at the following line number if the user forgets to set decomposition CalculateFacesSet to True.
https://github.com/OpenCMISS/iron/blob/devel/src/data_projection_routines.f90#L1111

Need to check if all projection types (ALL_ELEMENTS', 'BOUNDARY_FACES', 'BOUNDARY_LINES') use face calculations. If so, then an error should be raised if the user attempts to project data without setting decomposition CalculateFacesSet to True.

Currently the way to set candidate faces for data projection (DataProjection_ProjectionCandidatesSet) is by inputting an array of element numbers and another array of their corresponding local element faces. These candidate faces are set for all data points. There are several things to change about this:

1. We should have the option to be able to set different candidate faces for each data point.
2. The local face numbering order is currently -xi1, -xi2, -xi3, +xi1, +xi2, +xi3. It would make more sense to use -xi1, +xi1, -xi2, +xi2, -xi3, +xi3 instead.
3. Perhaps instead of specifying face numbers, which are heavily convention-dependent and therefore difficult for users who don't remember the convention, we could allow face specification using xi directions instead.

Also, it would be worthwhile to keep the face specification method consistent across different subroutines.

Some subroutines have a fairly large number of association and allocation pointer checks that are implemented in a very deeply nested IF - ELSE IF structure. This can be hard to hard for developers and possibly difficult for compilers to optimize (when dealing with a high level of nesting).

A possible alternative is to conduct the checks in a series of sequential IF-THEN statements at the beginning of the subroutine. For instance:

if ( .not.ASSOCIATED(DOMAIN) ) call FlagError( "Domain is not associated.", ERR,ERROR, *998 )
    if ( .not.ASSOCIATED(DOMAIN%MAPPINGS) ) call FlagError( "Domain mappings is not associated.", ERR,ERROR, *998 )
    if ( .not.ASSOCIATED(DOMAIN%MAPPINGS%ELEMENTS) ) call FlagError("Domain mappings elements is not associated.",ERR,ERROR,*998)
    if ( .not.ASSOCIATED(DOMAIN%DECOMPOSITION) ) call FlagError( "Domain decomposition is not associated.", ERR,ERROR, *998 )
if ( .not.ASSOCIATED(DOMAIN%MESH) ) call FlagError( "Domain mesh is not associated.", ERR,ERROR, *998) 

Not an essential issue but maybe a coding style we could adopt as code is added and modified.

There are two issues: backslashes in path are interpreted as escape characters, and on Windows dlls are expected to be in a \bin directory, not \lib.

This will be used with a meshfree implementation that Chun Goh is currently working on for his PhD.
Implementation notes described here:
https://docs.google.com/a/aucklanduni.ac.nz/document/d/1Rk7oMn2b-gGm_zWRUBtIRdC5HTeBCxRFIKq2TipcMfc/edit?usp=sharing

When doing data projection evaluations (DataProjection_DataPointsProjectionEvaluate), the DataProjection_Newton<Elements/Faces/Lines>Evaluate subroutines are all currently evaluating distance vectors by doing:

DISTANCE_VECTOR(1:REGION_DIMENSIONS)=POINT_VALUES-INTERPOLATED_POINT%VALUES(:,NO_PART_DERIV)

This is wrong because it means that the distance vectors are pointing from the projected coordinates to the data point coordinates, which is opposite to how it should be. While this doesn't affect the value of the distance calculation, it is counter to the concept of projection and it affects any visualisations of the projection vectors if we choose to output them.

Changing the order of subtraction would correct the problem however we must be careful to check whether the change in sign of the distance vector is going to affect minimisation maths associated with the projection functions.

Some of the tests for the CMake system seem to be been added in order to check all the outputs. Some of these tests have never worked and should probably not be tested. This is an issue for tests that should be removed.

• Darcy Exp and Darcy SinCos: These tests fail with a floating point exception. The reason for the floating point exception is that a variable, PERM_OVER_VIS, is NaN as it has never been initialised. I'm not sure what is going on with the code as it seems like it should never work. I'm not sure how to initialise the variable. We could, perhaps, contact Adam for comment but apart from that I'd suggest removing these tests.

At the moment the boundary conditions object is derived from the solver equations object. Now boundary conditions are technically applied to DOFs, that is a particular field DOF is flagged in some way e.g., fixed at a certain value. These fixed DOFs in turn affect the solution matrices i.e., rows and columns are eliminated from the equations matrices. This is why the BCs are set as part of the solver equations. However, for complicated workflows this is not so straight forward. I may have multiple solvers in a control loop. These solvers may each have a solver equations. These solver equations, however, may involve the same field variables. What if a BC was set on a DOF in one set of solver equations but not in the other? What if one BC value was set for a DOF in one equations set but a different BC value on the same DOF in the next equations set? Maybe the BC values should not be set as part of the fixed,free API and just take the value from current field variable value. Now it may be that we are doing some sort of iteration between different solvers and so we may want to fix and free the same DOF in different solves? Thus we do need a way of fixing and freeing DOFs on a solver equations level. I don't think it makes sense to have different DOF values??? There are also constraint conditions that we need to allow for e.g., linear constraints between DOFs. This issue is set up to discuss changes to the BC setting API.

There should be a CMake target "docs" that creates the documentation for the current Iron source.

Add in the ability to solve optimisation problems. Requires:

• Add in optimisation solver
• Add in interfaces to TAO solvers in PETSc
• Change equations to allow for three types of equations: scalar equations, vector equations (what we have now), functional equations (i.e., CellML type equations).
• Change equations to allow for different types of equalities. Equations are now = 0 but we will need > 0, < 0, >= 0 and <= 0 to allow for inequality constraints.
• For the scalar equations the mapping will be from field variable vectors to equations elements (c.f., dynamic, linear and non-linear elements in vector equations). Some ways to derive a scalar from the field variable vector include norms, vector dot products i.e., x^T.y, quadratic forms i.e., x^T.A.y etc. Will probably need to allow for all of these in addition to simple scalars and maybe a single scalar functional form (CellML)?
• Change solver equations to add in objective function (add a scalar equations set) and constraints (add in the vector equations sets to allow for PDE constrained optimisation, linear programs etc.). Possibly change SolverEquations to SolverSystemEquations and CellMLEquations to SolverFunctionalEquations etc.
• Might need to think that we have materials fields in the scalar equations set equations as well as dependent fields e.g., for a linear program the objective is c^T.x
• Add in Bound conditions on variables c.f., boundary and interface conditions. Will need to deal with equality and inequality bounds.
• Think about how gradients and Hessians are associated with the objective, constraint Jacobians etc.

From @chrispbradley on February 25, 2016 0:5

Not sure if this is a manage issue (i.e., cmake) or an iron issue but the iron visual studio project does not have the include files used in iron visible from visual studio (i.e., we should have an include folder as well as a src folder). Not sure if this means that the iron project does not depend on the include files or not? If not then it probably should as a change to an include file should trigger a re-build.

Copied from original issue: OpenCMISS/manage#60

It would be nice to have support for embedded meshes. A mesh should be able to store fields of another, finer mesh, that shares the geometric location. The positions of the coarse and fine nodes should be coupled, such that if the parent mesh deforms, also the positions of the sub mesh nodes are updated.
Also the decomposition could be the same for the coarse and fine elements.

Currently, a data point projection is created with a mesh. The geometric field for the projection is supplied when the projection is evaluated. The data point projection should be created with the field for the projection. The decomposition/mesh can then be obtained from this field. The data projection evaluate should take the parameter set of the field in order to evaluate the projection.

I've implemented this solver. I think it works but i need to test it properly, still.

also I implemented a routine to tell the forward and improved Euler solvers how many steps they should use instead of which step size. If one does not use the routine, everything will still work as before.

initially when I planned OpenCMISS I decided to, unlike old cmiss, not have the user access lines and faces via a global line or face number. Instead relative references via an element were envisaged e.g., give an element and then a relative reference to identify the desired line or face. Now, some of these references have crept in via local line and face numbers e.g., for a standard quad element there are 6 local faces and 12 local lines. These local references are, however, arbitrary e.g., what is local line 5? Instead the local lines and faces should be referenced via the xi directions normal to the desired line or face. For faces this is straight forward. It would be the index of the xi direction (outward) normal to the face. For quad elements this would be -3,-2,-1,+1,+2,+3 for the 6 faces. For tet elements the direction is given by the local node at the vertex opposite the face. For lines we need two normal xi directions to identify the two faces for which the line is the intersection of. Thus for a tri-linear quad element the xi 1 local line from local node 1 to 2 would be identified by the -3 and -2 xi direction indexes. Than plan would be to change all the APIs such that if the user needed to specify a local/global line or face they would use the element and xi normal indices. Any comments?

Currently no checks on the required python version are being done. this could also be an issue for the opencmiss-wide manage build system?!

It may be good to specify field variable attributes as part of the equations mapping. At the moment attributes like time dependency and linearity are deduced as part of the attributes of the equations that the variables are in. However, there are other variables for which this is not possible. For example, some equations sets have "auxiliary" field variables that are not directly included in the equations. At example might be stress and strain U1/U2 variables for elasticity. It would be good to be able to flag these variables as auxiliary and, say, time dependent so that they can be included in operations for other time dependent field variables e.g., copying to the previous value parameter set, writing them out as part of a history file.

The InterfaceMatrix_TimeDependenceTypeGet/Set routines are not in line with standard method naming conventions. There first argument is an interface condition type and not an interface matrix type. This causes a problem for the binding generation scripts which take their cues from the method name. As a result these routines are not available in python. They should be moved to the interfaceEquations methods and take the interfaceEquations object as their first parameter.

Issue for the discussion of what Fortran 2008 features can be used in OpenCMISS. In general, the use of language features depends on performance implications and the support of that feature by the various Fortran compilers. A useful summary of Fortran language features that have been implemented by the various compilers can be found here:

http://fortranwiki.org/fortran/show/Fortran+2008+status

OpenCMISS can currently interpolate strain at a given xi using Elasticity_StrainInterpolateXi.
However, this routines need to be generalised to output any relevant tensor (depends on the the equation set used).

See https://tracker.physiomeproject.org/show_bug.cgi?id=3126 for previous discussions on outputting stress and strain fields.