Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add new API for getting Jacobian information about model. #41

Merged
merged 8 commits into from
Feb 13, 2025
40 changes: 35 additions & 5 deletions src/scaffoldfitter/fitter.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,9 @@
from cmlibs.maths.vectorops import add, mult, sub
from cmlibs.utils.zinc.field import assignFieldParameters, createFieldFiniteElementClone, getGroupList, \
findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName
from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName, getMaximumNodeIdentifier
from cmlibs.utils.zinc.finiteelement import evaluate_field_nodeset_range, findNodeWithName
from cmlibs.utils.zinc.general import ChangeManager
from cmlibs.utils.zinc.mesh import calculate_jacobian, report_on_lowest_value
from cmlibs.utils.zinc.region import write_to_buffer, read_from_buffer
from cmlibs.zinc.context import Context
from cmlibs.zinc.element import Elementbasis, Elementfieldtemplate
Expand Down Expand Up @@ -186,7 +187,7 @@ def getInheritFitterStep(self, refFitterStep: FitterStep):
"""
refType = type(refFitterStep)
for index in range(self._fitterSteps.index(refFitterStep) - 1, -1, -1):
if type(self._fitterSteps[index]) == refType:
if type(self._fitterSteps[index]) is refType:
return self._fitterSteps[index]
return None

Expand Down Expand Up @@ -841,7 +842,7 @@ def getDataRMSAndMaximumProjectionErrorForGroup(self, groupName):
with ChangeManager(self._fieldmodule):
group = self._fieldmodule.findFieldByName(groupName).castGroup()
if group.isValid():
calculation_field = self._fieldmodule.createFieldAnd(group, self._activeDataGroupField)
calculation_field = self._fieldmodule.createFieldAnd(group, self._activeDataGroupField)
nodeset = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS)
temp_dataset_group = self._fieldmodule.createFieldGroup().createNodesetGroup(nodeset)
temp_dataset_group.addNodesConditional(calculation_field)
Expand All @@ -852,6 +853,37 @@ def getDataRMSAndMaximumProjectionErrorForGroup(self, groupName):

return None, None

def getModelWorstElementJacobianInfo(self, mesh_group=None):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

getLowestElementJacobian(self, mesh_group=None)
Add comment that a value <= 0.0 is bad, either zero or negative volume / left-handed coordinates.

"""
Get the information on the 3D element with the worst jacobian value (most negative).
Optional mesh group parameter allows the user to make the calculation over a different group
from the whole mesh.
:param mesh_group: Optional parameter to specify a particular mesh group to make the calculation over.
:return: Element identifier, minimum jacobian value. Values are -1, inf if there is no data or bad fields.
"""
with ChangeManager(self._fieldmodule):
jacobian = calculate_jacobian(self._modelCoordinatesField)
report = report_on_lowest_value(jacobian, mesh_group if mesh_group else None)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Add following line within context manager:
del jacobian
so no notification about field ever existing.
Note I'd also prefer mesh_3d to be discovered here and a mesh or mesh_group passed to the function creating the jacobian field.


return report

def getModelWorstElementJacobianInfoForGroup(self, group_name):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

getLowestElementJacobianForGroupName

"""
Get the information on the 3D element with the worst jacobian
value (most negative) of the 3D mesh group with the given name.
If the group_name is not a valid group name then None, None is returned.

:param group_name: Name of group to make calculation over.
:return: Element identifier, minimum jacobian value.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to note returns -1, inf if failed to evaluate, but this is at odds with returning None, None below.

"""
with ChangeManager(self._fieldmodule):
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ChangeManager is redundant here.

group = self._fieldmodule.findFieldByName(group_name).castGroup()
if group.isValid():
result = self.getMeshWorstElementJacobianInfo(group)
return result

return None, None

def getMarkerDataFields(self):
"""
Only call if markerGroup exists.
Expand Down Expand Up @@ -1096,7 +1128,6 @@ def _discoverModelFitGroup(self):
Discover modelFitGroup from set name.
"""
self._modelFitGroup = None
field = None
if self._modelFitGroupName:
field = self._fieldmodule.findFieldByName(self._modelFitGroupName)
if field.isValid():
Expand Down Expand Up @@ -1161,7 +1192,6 @@ def _discoverFlattenGroup(self):
Discover flattenGroup from set name.
"""
self._flattenGroup = None
field = None
if self._flattenGroupName:
field = self._fieldmodule.findFieldByName(self._flattenGroupName)
if field.isValid():
Expand Down