diff --git a/README.rst b/README.rst new file mode 100644 index 0000000..7ea7e24 --- /dev/null +++ b/README.rst @@ -0,0 +1,13 @@ +Scaffold Fitter +=============== + +Scaffold/model fitter using OpenCMISS-Zinc. + +For the interim, install with the following command:: + + pip install -e . + +Requires the following Python libraries to be installed: + +- opencmiss.zinc (manual install from opencmiss.org) +- github.com/OpenCMISS-Bindings/opencmiss.utils diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..bee40e4 --- /dev/null +++ b/setup.py @@ -0,0 +1,45 @@ +from setuptools import setup, find_packages +import os +import io + +SETUP_DIR = os.path.dirname(os.path.abspath(__file__)) + +# List all of your Python package dependencies in the +# requirements.txt file + +def readfile(filename, split=False): + with io.open(filename, encoding="utf-8") as stream: + if split: + return stream.read().split("\n") + return stream.read() + +readme = readfile("README.rst", split=True)[3:] # skip title +# For requirements not hosted on PyPi place listings +# into the 'requirements.txt' file. +requires = [ + # minimal requirements listing + "opencmiss.utils >= 0.2", + "opencmiss.zinc" # not yet on pypi - need manual install from opencmiss.org +] +source_license = readfile("LICENSE") + +setup( + name="scaffoldfitter", + version="0.1.1", + description="Scaffold/model geometric fitting library using OpenCMISS-Zinc.", + long_description="\n".join(readme) + source_license, + classifiers=[ + "Development Status :: 3 - Alpha", + "Programming Language :: Python", + "Topic :: Scientific/Engineering :: Medical Science Apps." + ], + author="Auckland Bioengineering Institute", + author_email="r.christie@auckland.ac.nz", + url="https://github.com/ABI-Software/scaffoldfitter", + license="Apache Software License", + packages=find_packages("src"), + package_dir={"": "src"}, + include_package_data=True, + zip_safe=False, + install_requires=requires, + ) diff --git a/src/scaffoldfitter/__init__.py b/src/scaffoldfitter/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/src/scaffoldfitter/fitter.py b/src/scaffoldfitter/fitter.py new file mode 100644 index 0000000..290636e --- /dev/null +++ b/src/scaffoldfitter/fitter.py @@ -0,0 +1,717 @@ +""" +Main class for fitting scaffolds. +""" + +import json +from opencmiss.utils.zinc.field import assignFieldParameters, createFieldFiniteElementClone, getGroupList, getManagedFieldNames, \ + findOrCreateFieldFiniteElement, findOrCreateFieldStoredMeshLocation, getUniqueFieldName, orphanFieldByName +from opencmiss.utils.zinc.finiteelement import evaluateFieldNodesetMean, findNodeWithName, getMaximumNodeIdentifier +from opencmiss.utils.zinc.general import ChangeManager +from opencmiss.zinc.context import Context +from opencmiss.zinc.field import Field, FieldFindMeshLocation, FieldGroup +from opencmiss.zinc.result import RESULT_OK, RESULT_WARNING_PART_DONE + + +class Fitter: + + def __init__(self, zincModelFileName : str, zincDataFileName : str): + """ + :param zincModelFileName: Name of zinc file supplying model to fit. + :param zincDataFileName: Name of zinc filed supplying data to fit to. + """ + self._zincModelFileName = zincModelFileName + self._zincDataFileName = zincDataFileName + self._context = Context("Scaffoldfitter") + self._region = None + self._rawDataRegion = None + self._fieldmodule = None + self._modelCoordinatesField = None + self._modelCoordinatesFieldName = None + self._modelReferenceCoordinatesField = None + self._dataCoordinatesField = None + self._dataCoordinatesFieldName = None + self._mesh = [] # [dimension - 1] + self._dataProjectionLocationFields = [ ] # [dimension - 1] + self._dataProjectionCoordinatesFields = [ ] # [dimension - 1] + self._dataProjectionDeltaFields = [ ] # [dimension - 1] + self._dataProjectionErrorFields = [ ] # [dimension - 1] + self._dataProjectionNodeGroupFields = [] # [dimension - 1] + self._dataProjectionNodesetGroups = [] # [dimension - 1] + self._dataProjectionDirectionField = None # for storing original projection direction unit vector + self._markerGroup = None + self._markerGroupName = None + self._markerNodeGroup = None + self._markerLocationField = None + self._markerNameField = None + self._markerCoordinatesField = None + self._markerDataGroup = None + self._markerDataCoordinatesField = None + self._markerDataNameField = None + self._markerDataLocationField = None + self._markerDataLocationCoordinatesField = None + self._markerDataDeltaField = None + self._markerDataLocationGroupField = None + self._markerDataLocationGroup = None + self._diagnosticLevel = 0 + self._fitterSteps = [] + + def decodeSettingsJSON(self, s : str, decoder): + """ + Define Fitter from JSON serialisation output by encodeSettingsJSON. + :param s: String of JSON encoded Fitter settings. + :param decoder: decodeJSONFitterSteps(fitter, dct) for decodings FitterSteps. + """ + self._fitterSteps.clear() + dct = json.loads(s, object_hook=lambda dct: decoder(self, dct)) + self._modelCoordinatesFieldName = dct["modelCoordinatesField"] + self._dataCoordinatesFieldName = dct["dataCoordinatesField"] + self._markerGroupName = dct["markerGroup"] + self._diagnosticLevel = dct["diagnosticLevel"] + # self._fitterSteps will already be populated by decoder so don't need + #self._fitterSteps = dct["fitterSteps"] + + def encodeSettingsJSON(self) -> str: + """ + :return: String JSON encoding of Fitter settings. + """ + dct = { + "modelCoordinatesField" : self._modelCoordinatesFieldName, + "dataCoordinatesField" : self._dataCoordinatesFieldName, + "markerGroup" : self._markerGroupName, + "diagnosticLevel" : self._diagnosticLevel, + "fitterSteps" : [ fitterStep.encodeSettingsJSONDict() for fitterStep in self._fitterSteps ] + } + return json.dumps(dct, sort_keys=False, indent=4) + + def load(self): + """ + Read model and data and define fit fields and data. + Can call again to reset fit, after parameters have changed. + """ + self._region = self._context.createRegion() + self._fieldmodule = self._region.getFieldmodule() + self._rawDataRegion = self._region.createChild("raw_data") + self._loadModel() + self._loadData() + self._defineDataProjectionFields() + self.calculateDataProjections() + + def _loadModel(self): + result = self._region.readFile(self._zincModelFileName) + assert result == RESULT_OK, "Failed to load model file" + str(self._zincModelFileName) + self._mesh = [ self._fieldmodule.findMeshByDimension(d + 1) for d in range(3) ] + self._discoverModelCoordinatesField() + + def _loadData(self): + """ + Load zinc data file into self._rawDataRegion. + Transfer data points (and converted nodes) into self._region. + """ + result = self._rawDataRegion.readFile(self._zincDataFileName) + assert result == RESULT_OK, "Failed to load data file " + str(self._zincDataFileName) + # if there both nodes and datapoints, offset datapoint identifiers to ensure different + fieldmodule = self._rawDataRegion.getFieldmodule() + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + if nodes.getSize() > 0: + datapoints = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + if datapoints.getSize() > 0: + maximumDatapointIdentifier = max(0, getMaximumNodeIdentifier(datapoints)) + maximumNodeIdentifier = max(0, getMaximumNodeIdentifier(nodes)) + # this assumes identifiers are in low ranges and can be improved if there is a problem: + identifierOffset = 100000 + while (maximumDatapointIdentifier > identifierOffset) or (maximumNodeIdentifier > identifierOffset): + assert identifierOffset < 1000000000, "Invalid node and datapoint identifier ranges" + identifierOffset *= 10 + with ChangeManager(fieldmodule): + while True: + # logic relies on datapoints being in identifier order + datapoint = datapoints.createNodeiterator().next() + identifier = datapoint.getIdentifier() + if identifier >= identifierOffset: + break; + result = datapoint.setIdentifier(identifier + identifierOffset) + assert result == RESULT_OK, "Failed to offset datapoint identifier" + # transfer nodes as datapoints to self._region + sir = self._rawDataRegion.createStreaminformationRegion() + srm = sir.createStreamresourceMemory() + sir.setResourceDomainTypes(srm, Field.DOMAIN_TYPE_NODES) + self._rawDataRegion.write(sir) + result, buffer = srm.getBuffer() + assert result == RESULT_OK, "Failed to write nodes" + buffer = buffer.replace(bytes("!#nodeset nodes", "utf-8"), bytes("!#nodeset datapoints", "utf-8")) + sir = self._region.createStreaminformationRegion() + srm = sir.createStreamresourceMemoryBuffer(buffer) + result = self._region.read(sir) + assert result == RESULT_OK, "Failed to load nodes as datapoints" + # transfer datapoints to self._region + sir = self._rawDataRegion.createStreaminformationRegion() + srm = sir.createStreamresourceMemory() + sir.setResourceDomainTypes(srm, Field.DOMAIN_TYPE_DATAPOINTS) + self._rawDataRegion.write(sir) + result, buffer = srm.getBuffer() + assert result == RESULT_OK, "Failed to write datapoints" + sir = self._region.createStreaminformationRegion() + srm = sir.createStreamresourceMemoryBuffer(buffer) + result = self._region.read(sir) + assert result == RESULT_OK, "Failed to load datapoints" + self._discoverDataCoordinatesField() + self._discoverMarkerGroup() + + def getDataCoordinatesField(self): + return self._dataCoordinatesField + + def setDataCoordinatesField(self, dataCoordinatesField : Field): + finiteElementField = dataCoordinatesField.castFiniteElement() + assert finiteElementField.isValid() and (finiteElementField.getNumberOfComponents() == 3) + self._dataCoordinatesFieldName = dataCoordinatesField.getName() + self._dataCoordinatesField = finiteElementField + + def setDataCoordinatesFieldByName(self, dataCoordinatesFieldName): + self.setDataCoordinatesField(self._fieldmodule.findFieldByName(dataCoordinatesFieldName)) + + def _discoverDataCoordinatesField(self): + """ + Choose default dataCoordinates field. + """ + self._dataCoordinatesField = None + field = None + if self._dataCoordinatesFieldName: + field = self._fieldmodule.findFieldByName(self._dataCoordinatesFieldName) + else: + datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + datapoint = datapoints.createNodeiterator().next() + if datapoint.isValid(): + fieldcache = self._fieldmodule.createFieldcache() + fieldcache.setNode(datapoint) + fielditer = self._fieldmodule.createFielditerator() + field = fielditer.next() + while field.isValid(): + if field.isTypeCoordinate() and (field.getNumberOfComponents() == 3) and (field.castFiniteElement().isValid()): + if field.isDefinedAtLocation(fieldcache): + break; + field = fielditer.next() + else: + field = None + self.setDataCoordinatesField(field) + + def getMarkerGroup(self): + return self._markerGroup + + def setMarkerGroup(self, markerGroup : Field): + self._markerGroup = None + self._markerGroupName = None + self._markerNodeGroup = None + self._markerLocationField = None + self._markerCoordinatesField = None + self._markerNameField = None + self._markerDataGroup = None + self._markerDataCoordinatesField = None + self._markerDataNameField = None + self._markerDataLocationField = None + self._markerDataLocationCoordinatesField = None + self._markerDataDeltaField = None + self._markerDataLocationGroupField = None + self._markerDataLocationGroup = None + if not markerGroup: + return + fieldGroup = markerGroup.castGroup() + assert fieldGroup.isValid() + self._markerGroup = fieldGroup + self._markerGroupName = markerGroup.getName() + nodes = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + self._markerNodeGroup = self._markerGroup.getFieldNodeGroup(nodes).getNodesetGroup() + if self._markerNodeGroup.isValid(): + node = self._markerNodeGroup.createNodeiterator().next() + if node.isValid(): + fieldcache = self._fieldmodule.createFieldcache() + fieldcache.setNode(node) + fielditer = self._fieldmodule.createFielditerator() + field = fielditer.next() + while field.isValid(): + if field.isDefinedAtLocation(fieldcache): + if (not self._markerLocationField) and field.castStoredMeshLocation().isValid(): + self._markerLocationField = field + elif (not self._markerNameField) and (field.getValueType() == Field.VALUE_TYPE_STRING): + self._markerNameField = field + field = fielditer.next() + self._updateMarkerCoordinatesField() + else: + self._markerNodeGroup = None + datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + self._markerDataGroup = self._markerGroup.getFieldNodeGroup(datapoints).getNodesetGroup() + if self._markerDataGroup.isValid(): + datapoint = self._markerDataGroup.createNodeiterator().next() + if datapoint.isValid(): + fieldcache = self._fieldmodule.createFieldcache() + fieldcache.setNode(datapoint) + fielditer = self._fieldmodule.createFielditerator() + field = fielditer.next() + while field.isValid(): + if field.isDefinedAtLocation(fieldcache): + if (not self._markerDataCoordinatesField) and field.isTypeCoordinate() and \ + (field.getNumberOfComponents() == 3) and (field.castFiniteElement().isValid()): + self._markerDataCoordinatesField = field + elif (not self._markerDataNameField) and (field.getValueType() == Field.VALUE_TYPE_STRING): + self._markerDataNameField = field + field = fielditer.next() + else: + self._markerDataGroup = None + self._calculateMarkerDataLocations() + + def setMarkerGroupByName(self, markerGroupName): + self.setMarkerGroup(self._fieldmodule.findFieldByName(markerGroupName)) + + def getMarkerDataFields(self): + """ + Only call if markerGroup exists. + :return: markerDataGroup, markerDataCoordinates, markerDataName + """ + return self._markerDataGroup, self._markerDataCoordinatesField, self._markerDataNameField + + def getMarkerDataLocationFields(self): + """ + Get fields giving marker location coordinates and delta on the data points (copied from nodes). + Only call if markerGroup exists. + :return: markerDataLocation, markerDataLocationCoordinates, markerDataDelta + """ + return self._markerDataLocationField, self._markerDataLocationCoordinatesField, self._markerDataDeltaField + + def getMarkerModelFields(self): + """ + Only call if markerGroup exists. + :return: markerNodeGroup, markerLocation, markerCoordinates, markerName + """ + return self._markerNodeGroup, self._markerLocationField, self._markerCoordinatesField, self._markerNameField + + def _calculateMarkerDataLocations(self): + """ + Called when markerGroup exists. + Find matching marker mesh locations for marker data points. + Only finds matching location where there is one datapoint and one node + for each name in marker group. + Defines datapoint group self._markerDataLocationGroup to contain those with locations. + """ + self._markerDataLocationField = None + self._markerDataLocationCoordinatesField = None + self._markerDataDeltaField = None + self._markerDataLocationGroupField = None + self._markerDataLocationGroup = None + if not (self._markerDataGroup and self._markerDataNameField and self._markerNodeGroup and self._markerLocationField and self._markerNameField): + return + + markerPrefix = self._markerGroup.getName() + # assume marker locations are in highest dimension mesh + mesh = self.getHighestDimensionMesh() + datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + meshDimension = mesh.getDimension() + fieldcache = self._fieldmodule.createFieldcache() + with ChangeManager(self._fieldmodule): + self._markerDataLocationField = findOrCreateFieldStoredMeshLocation(self._fieldmodule, mesh, name=markerPrefix + "_data_location_" + mesh.getName(), managed=False) + self._markerDataLocationGroupField = self._fieldmodule.createFieldNodeGroup(datapoints) + self._markerDataLocationGroupField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_data_location_group")) + self._markerDataLocationGroup = self._markerDataLocationGroupField.getNodesetGroup() + self._updateMarkerDataLocationCoordinatesField() + nodetemplate = self._markerDataGroup.createNodetemplate() + nodetemplate.defineField(self._markerDataLocationField) + datapointIter = self._markerDataGroup.createNodeiterator() + datapoint = datapointIter.next() + while datapoint.isValid(): + fieldcache.setNode(datapoint) + name = self._markerDataNameField.evaluateString(fieldcache) + # if this is the only datapoint with name: + if name and findNodeWithName(self._markerDataGroup, self._markerDataNameField, name): + node = findNodeWithName(self._markerNodeGroup, self._markerNameField, name) + if node: + fieldcache.setNode(node) + element, xi = self._markerLocationField.evaluateMeshLocation(fieldcache, meshDimension) + if element.isValid(): + datapoint.merge(nodetemplate) + fieldcache.setNode(datapoint) + self._markerDataLocationField.assignMeshLocation(fieldcache, element, xi) + self._markerDataLocationGroup.addNode(datapoint) + datapoint = datapointIter.next() + # Warn about datapoints without a location in model + markerDataGroupSize = self._markerDataGroup.getSize() + markerDataLocationGroupSize = self._markerDataLocationGroup.getSize() + markerNodeGroupSize = self._markerNodeGroup.getSize() + if self.getDiagnosticLevel() > 0: + if markerDataLocationGroupSize < markerDataGroupSize: + print("Warning: Only " + str(markerDataLocationGroupSize) + " of " + str(markerDataGroupSize) + " marker data points have model locations") + if markerDataLocationGroupSize < markerNodeGroupSize: + print("Warning: Only " + str(markerDataLocationGroupSize) + " of " + str(markerNodeGroupSize) + " marker model locations used") + + def _discoverMarkerGroup(self): + self._markerGroup = None + self._markerNodeGroup = None + self._markerLocationField = None + self._markerNameField = None + self._markerCoordinatesField = None + markerGroup = self._fieldmodule.findFieldByName(self._markerGroupName if self._markerGroupName else "marker").castGroup() + if not markerGroup.isValid(): + markerGroup = None + self.setMarkerGroup(markerGroup) + + def _updateMarkerCoordinatesField(self): + if self._modelCoordinatesField and self._markerLocationField: + with ChangeManager(self._fieldmodule): + markerPrefix = self._markerGroup.getName() + self._markerCoordinatesField = self._fieldmodule.createFieldEmbedded(self._modelCoordinatesField, self._markerLocationField) + self._markerCoordinatesField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_coordinates")) + else: + self._markerCoordinatesField = None + + def _updateMarkerDataLocationCoordinatesField(self): + if self._modelCoordinatesField and self._markerDataLocationField: + with ChangeManager(self._fieldmodule): + markerPrefix = self._markerGroup.getName() + self._markerDataLocationCoordinatesField = self._fieldmodule.createFieldEmbedded(self._modelCoordinatesField, self._markerDataLocationField) + self._markerDataLocationCoordinatesField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_data_location_coordinates")) + self._markerDataDeltaField = self._markerDataLocationCoordinatesField - self._markerDataCoordinatesField + self._markerDataDeltaField.setName(getUniqueFieldName(self._fieldmodule, markerPrefix + "_data_delta")) + else: + self._markerDataLocationCoordinatesField = None + + def getModelCoordinatesField(self): + return self._modelCoordinatesField + + def getModelReferenceCoordinatesField(self): + return self._modelReferenceCoordinatesField + + def setModelCoordinatesField(self, modelCoordinatesField : Field): + finiteElementField = modelCoordinatesField.castFiniteElement() + assert finiteElementField.isValid() and (finiteElementField.getNumberOfComponents() == 3) + self._modelCoordinatesField = finiteElementField + self._modelCoordinatesFieldName = modelCoordinatesField.getName() + modelReferenceCoordinatesFieldName = "reference_" + self._modelCoordinatesField.getName(); + orphanFieldByName(self._fieldmodule, modelReferenceCoordinatesFieldName); + self._modelReferenceCoordinatesField = createFieldFiniteElementClone(self._modelCoordinatesField, modelReferenceCoordinatesFieldName) + self._updateMarkerCoordinatesField() + self._updateMarkerDataLocationCoordinatesField() + + def setModelCoordinatesFieldByName(self, modelCoordinatesFieldName): + self.setModelCoordinatesField(self._fieldmodule.findFieldByName(modelCoordinatesFieldName)) + + def _discoverModelCoordinatesField(self): + """ + Choose default modelCoordinates field. + """ + self._modelCoordinatesField = None + self._modelReferenceCoordinatesField = None + field = None + if self._modelCoordinatesFieldName: + field = self._fieldmodule.findFieldByName(self._modelCoordinatesFieldName) + else: + mesh = self.getHighestDimensionMesh() + element = mesh.createElementiterator().next() + if element.isValid(): + fieldcache = self._fieldmodule.createFieldcache() + fieldcache.setElement(element) + fielditer = self._fieldmodule.createFielditerator() + field = fielditer.next() + while field.isValid(): + if field.isTypeCoordinate() and (field.getNumberOfComponents() == 3) and (field.castFiniteElement().isValid()): + if field.isDefinedAtLocation(fieldcache): + break; + field = fielditer.next() + else: + field = None + if field: + self.setModelCoordinatesField(field) + + def runConfig(self): + """ + Complete configuration setup. + """ + self._hasRunConfig = True + + def _defineDataProjectionFields(self): + self._dataProjectionLocationFields = [] + self._dataProjectionCoordinatesFields = [] + self._dataProjectionDeltaFields = [] + self._dataProjectionErrorFields = [] + self._dataProjectionNodeGroupFields = [] + self._dataProjectionNodesetGroups = [] + with ChangeManager(self._fieldmodule): + datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + for d in range(2): + mesh = self._mesh[d] + dataProjectionLocation = findOrCreateFieldStoredMeshLocation(self._fieldmodule, mesh, name="data_projection_location_" + mesh.getName(), managed=False) + self._dataProjectionLocationFields.append(dataProjectionLocation) + dataProjectionCoordinates = self._fieldmodule.createFieldEmbedded(self._modelCoordinatesField, dataProjectionLocation) + dataProjectionCoordinates.setName(getUniqueFieldName(self._fieldmodule, "data_projection_coordinates_" + mesh.getName())) + self._dataProjectionCoordinatesFields.append(dataProjectionCoordinates) + dataProjectionDelta = dataProjectionCoordinates - self._dataCoordinatesField + dataProjectionDelta.setName(getUniqueFieldName(self._fieldmodule, "data_projection_delta_" + mesh.getName())) + self._dataProjectionDeltaFields.append(dataProjectionDelta) + dataProjectionError = self._fieldmodule.createFieldMagnitude(dataProjectionDelta) + dataProjectionError.setName(getUniqueFieldName(self._fieldmodule, "data_projection_error_" + mesh.getName())) + self._dataProjectionErrorFields.append(dataProjectionError) + field = self._fieldmodule.createFieldNodeGroup(datapoints) + field.setName(getUniqueFieldName(self._fieldmodule, "data_projection_group_" + mesh.getName())) + self._dataProjectionNodeGroupFields.append(field) + self._dataProjectionNodesetGroups.append(field.getNodesetGroup()) + self._dataProjectionDirectionField = findOrCreateFieldFiniteElement(self._fieldmodule, "data_projection_direction", + components_count = 3, component_names = [ "x", "y", "z" ]) + + def calculateDataProjections(self): + """ + Find projections of datapoints' coordinates onto model coordinates, + by groups i.e. from datapoints group onto matching 2-D or 1-D mesh group. + Calculate and store projection direction unit vector. + """ + assert self._dataCoordinatesField and self._modelCoordinatesField + with ChangeManager(self._fieldmodule): + findMeshLocation = None + datapoints = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS) + fieldcache = self._fieldmodule.createFieldcache() + for d in range(2): + self._dataProjectionNodesetGroups[d].removeAllNodes() + groups = getGroupList(self._fieldmodule) + for group in groups: + groupName = group.getName() + dataGroup = group.getFieldNodeGroup(datapoints).getNodesetGroup() + if not dataGroup.isValid(): + continue + for dimension in range(2, 0, -1): + meshGroup = group.getFieldElementGroup(self._mesh[dimension - 1]).getMeshGroup() + if meshGroup.isValid() and (meshGroup.getSize() > 0): + break + else: + if self.getDiagnosticLevel() > 0: + print("Fit Geometry: Warning: Cannot project data for group " + groupName + " as no matching mesh group") + continue + meshLocation = self._dataProjectionLocationFields[dimension - 1] + dataProjectionNodesetGroup = self._dataProjectionNodesetGroups[dimension - 1] + nodeIter = dataGroup.createNodeiterator() + node = nodeIter.next() + fieldcache.setNode(node) + if not self._dataCoordinatesField.isDefinedAtLocation(fieldcache): + if self.getDiagnosticLevel() > 0: + print("Fit Geometry: Warning: Cannot project data for group " + groupName + " as field " + self._dataCoordinatesField.getName() + " is not defined on data") + continue + if not meshLocation.isDefinedAtLocation(fieldcache): + # define meshLocation and on data Group: + nodetemplate = datapoints.createNodetemplate() + nodetemplate.defineField(meshLocation) + nodetemplate.defineField(self._dataProjectionDirectionField) + while node.isValid(): + result = node.merge(nodetemplate) + #print("node",node.getIdentifier(),"result",result) + node = nodeIter.next() + del nodetemplate + # restart iteration + nodeIter = dataGroup.createNodeiterator() + node = nodeIter.next() + findMeshLocation = self._fieldmodule.createFieldFindMeshLocation(self._dataCoordinatesField, self._modelCoordinatesField, meshGroup) + findMeshLocation.setSearchMode(FieldFindMeshLocation.SEARCH_MODE_NEAREST) + while node.isValid(): + fieldcache.setNode(node) + element, xi = findMeshLocation.evaluateMeshLocation(fieldcache, dimension) + if not element.isValid(): + print("Fit Geometry: Error finding data projection nearest mesh location for group " + groupName + ". Aborting group.") + break + result = meshLocation.assignMeshLocation(fieldcache, element, xi) + #print(result, "node", node.getIdentifier(), "element", element.getIdentifier(), "xi", xi) + #if result != RESULT_OK: + # mesh = meshLocation.getMesh() + # print("--> mesh", mesh.isValid(), mesh.getDimension(), findMeshLocation.getMesh().getDimension()) + # print("node", node.getIdentifier(), "is defined", meshLocation.isDefinedAtLocation(fieldcache)) + assert result == RESULT_OK, "Fit Geometry: Failed to assign data projection mesh location for group " + groupName + dataProjectionNodesetGroup.addNode(node) + node = nodeIter.next() + + # Store data projection directions + for dimension in range(1, 3): + nodesetGroup = self._dataProjectionNodesetGroups[dimension - 1] + if nodesetGroup.getSize() > 0: + fieldassignment = self._dataProjectionDirectionField.createFieldassignment( + self._fieldmodule.createFieldNormalise(self._dataProjectionDeltaFields[dimension - 1])) + fieldassignment.setNodeset(nodesetGroup) + result = fieldassignment.assign() + assert result in [ RESULT_OK, RESULT_WARNING_PART_DONE ], \ + "Fit Geometry: Failed to assign data projection directions for dimension " + str(dimension) + + if self.getDiagnosticLevel() > 0: + # Warn about unprojected points + unprojectedDatapoints = self._fieldmodule.createFieldNodeGroup(datapoints).getNodesetGroup() + unprojectedDatapoints.addNodesConditional(self._fieldmodule.createFieldIsDefined(self._dataCoordinatesField)) + for d in range(2): + unprojectedDatapoints.removeNodesConditional(self._dataProjectionNodeGroupFields[d]) + unprojectedCount = unprojectedDatapoints.getSize() + if unprojectedCount > 0: + print("Warning: " + str(unprojected) + " data points with data coordinates have not been projected") + del unprojectedDatapoints + + # remove temporary objects before ChangeManager exits + del findMeshLocation + del fieldcache + + def _addFitterStep(self, fitterStep): + self._fitterSteps.append(fitterStep) + + def _removeFitterStep(self, fitterStep): + self._fitterSteps.remove(fitterStep) + + def getNextFitterStep(self, refFitterStep): + """ + Return next fitter step after refFitterStep, or before if last, otherwise None. + """ + index = self._fitterSteps.index(refFitterStep) + 1 + if index >= len(self._fitterSteps): + index -= 2 + if index < 0: + return None + return self._fitterSteps[index] + + def getDataProjectionDirectionField(self): + return self._dataProjectionDirectionField + + def getDataProjectionNodeGroupField(self, dimension): + assert 1 <= dimension <= 2 + return self._dataProjectionNodeGroupFields[dimension - 1] + + def getDataProjectionNodesetGroup(self, dimension): + assert 1 <= dimension <= 2 + return self._dataProjectionNodesetGroups[dimension - 1] + + def getDataProjectionLocationField(self, dimension): + """ + :return: Field giving element:xi location for data points on mesh of dimension. + """ + assert 1 <= dimension <= 2 + return self._dataProjectionLocationFields[dimension - 1] + + def getDataProjectionCoordinatesField(self, dimension): + """ + :return: Field giving coordinates of projections of data points on mesh of dimension. + """ + assert 1 <= dimension <= 2 + return self._dataProjectionCoordinatesFields[dimension - 1] + + def getDataProjectionDeltaField(self, dimension): + """ + :return: Field giving delta coordinates (projection coordinates - data coordinates) + for data points on mesh of dimension. + """ + assert 1 <= dimension <= 2 + return self._dataProjectionDeltaFields[dimension - 1] + + def getDataProjectionErrorField(self, dimension): + """ + :return: Field giving magnitude of data point delta coordinates. + """ + assert 1 <= dimension <= 2 + return self._dataProjectionErrorFields[dimension - 1] + + def getMarkerDataLocationGroupField(self): + return self._markerDataLocationGroupField + + def getMarkerDataLocationNodesetGroup(self): + return self._markerDataLocationGroup + + def getMarkerDataLocationField(self): + return self._markerDataLocationField + + def getContext(self): + return self._context + + def getRegion(self): + return self._region + + def getFieldmodule(self): + return self._fieldmodule + + def getFitterSteps(self): + return self._fitterSteps + + def getMesh(self, dimension): + assert 1 <= dimension <= 3 + return self._mesh[dimension - 1] + + def getHighestDimensionMesh(self): + """ + :return: Highest dimension mesh with elements in it, or None if none. + """ + for d in range(2, -1, -1): + mesh = self._mesh[d] + if mesh.getSize() > 0: + return mesh + return None + + def evaluateNodeGroupMeanCoordinates(self, groupName, coordinatesFieldName, isData = False): + group = self._fieldmodule.findFieldByName(groupName).castGroup() + assert group.isValid() + nodeset = self._fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_DATAPOINTS if isData else Field.DOMAIN_TYPE_NODES) + nodesetGroup = group.getFieldNodeGroup(nodeset).getNodesetGroup() + assert nodesetGroup.isValid() + coordinates = self._fieldmodule.findFieldByName(coordinatesFieldName) + return evaluateFieldNodesetMean(coordinates, nodesetGroup) + + def getDiagnosticLevel(self): + return self._diagnosticLevel + + def setDiagnosticLevel(self, diagnosticLevel): + """ + :param diagnosticLevel: 0 = no diagnostic messages. 1 = Information and warning messages. 2 = Also optimisation reports. + """ + assert diagnosticLevel >= 0 + self._diagnosticLevel = diagnosticLevel + + def updateModelReferenceCoordinates(self): + assignFieldParameters(self._modelReferenceCoordinatesField, self._modelCoordinatesField) + + def writeModel(self, modelFileName): + """ + Write model nodes and elements excluding unmanaged fields to file. + """ + sir = self._region.createStreaminformationRegion() + srf = sir.createStreamresourceFile(modelFileName) + sir.setResourceFieldNames(srf, getManagedFieldNames(self._fieldmodule)) + sir.setResourceDomainTypes(srf, Field.DOMAIN_TYPE_NODES | Field.DOMAIN_TYPE_MESH1D | Field.DOMAIN_TYPE_MESH2D | Field.DOMAIN_TYPE_MESH3D) + result = self._region.write(sir) + assert result == RESULT_OK + + def writeData(self, fileName): + sir = self._region.createStreaminformationRegion() + sr = sir.createStreamresourceFile(fileName) + sir.setResourceDomainTypes(sr, Field.DOMAIN_TYPE_DATAPOINTS) + self._region.write(sir) + + +class FitterStep: + """ + Base class for fitter steps. + """ + + def __init__(self, fitter : Fitter): + """ + Construct and add to Fitter. + """ + self._fitter = fitter + fitter._addFitterStep(self) + self._hasRun = False + + def destroy(self): + """ + Remove from Fitter. + """ + self._fitter._removeFitterStep(self) + self._fitter = None + + def getFitter(self): + return self._fitter + + def hasRun(self): + return self._hasRun + + def setHasRun(self, hasRun): + self._hasRun = hasRun + + def getDiagnosticLevel(self): + return self._fitter.getDiagnosticLevel() + + def run(self): + """ + Override to perform action of derived FitStep + """ + pass diff --git a/src/scaffoldfitter/fitterstepalign.py b/src/scaffoldfitter/fitterstepalign.py new file mode 100644 index 0000000..b33078d --- /dev/null +++ b/src/scaffoldfitter/fitterstepalign.py @@ -0,0 +1,210 @@ +""" +Fit step for gross alignment and scale. +""" + +from opencmiss.utils.zinc.field import assignFieldParameters, createFieldsTransformations +from opencmiss.utils.zinc.finiteelement import getNodeNameCentres +from opencmiss.utils.zinc.general import ChangeManager +from opencmiss.zinc.field import Field +from opencmiss.zinc.optimisation import Optimisation +from opencmiss.zinc.result import RESULT_OK, RESULT_WARNING_PART_DONE +from scaffoldfitter.fitter import Fitter, FitterStep + +class FitterStepAlign(FitterStep): + + _jsonTypeId = "_FitterStepAlign" + + def __init__(self, fitter : Fitter): + super(FitterStepAlign, self).__init__(fitter) + self._alignMarkers = False + markerNodeGroup, markerLocation, markerCoordinates, markerName = fitter.getMarkerModelFields() + markerDataGroup, markerDataCoordinates, markerDataName = fitter.getMarkerDataFields() + if markerNodeGroup and markerLocation and markerCoordinates and markerName and \ + markerDataGroup and markerDataCoordinates and markerDataName: + self._alignMarkers = True + self._rotation = [ 0.0, 0.0, 0.0 ] + self._scale = 1.0 + self._translation = [ 0.0, 0.0, 0.0 ] + + @classmethod + def getJsonTypeId(cls): + return cls._jsonTypeId + + def decodeSettingsJSONDict(self, dct : dict): + """ + Decode definition of step from JSON dict. + """ + assert self._jsonTypeId in dct + self._alignMarkers = dct["alignMarkers"] + self._rotation = dct["rotation"] + self._scale = dct["scale"] + self._translation = dct["translation"] + + def encodeSettingsJSONDict(self) -> dict: + """ + Encode definition of step in dict. + :return: Settings in a dict ready for passing to json.dump. + """ + return { + self._jsonTypeId : True, + "alignMarkers" : self._alignMarkers, + "rotation" : self._rotation, + "scale" : self._scale, + "translation" : self._translation + } + + def isAlignMarkers(self): + return self._alignMarkers + + def setAlignMarkers(self, alignMarkers): + """ + :param alignMarkers: True to automatically align to markers, otherwise False. + """ + self._alignMarkers = alignMarkers + + def getRotation(self): + return self._rotation + + def setRotation(self, rotation): + """ + :param rotation: List of 3 euler angles in radians, order applied: + 0 = azimuth (about z) + 1 = elevation (about rotated y) + 2 = roll (about rotated x) + """ + assert len(rotation) == 3, "FitterStepAlign: Invalid rotation" + self._rotation = rotation + + def getScale(self): + return self._scale + + def setScale(self, scale): + """ + :param scale: Real scale. + """ + self._scale = scale + + def getTranslation(self): + return self._translation + + def setTranslation(self, translation): + """ + :param translation: [ x, y, z ]. + """ + assert len(translation) == 3, "FitterStepAlign: Invalid translation" + self._translation = translation + + def run(self): + """ + Perform align and scale. + """ + modelCoordinates = self._fitter.getModelCoordinatesField() + assert modelCoordinates, "Align: Missing model coordinates" + if self._alignMarkers: + self._doAlignMarkers() + fieldmodule = self._fitter._fieldmodule + with ChangeManager(fieldmodule): + # rotate, scale and translate model + modelCoordinatesTransformed = createFieldsTransformations( + modelCoordinates, self._rotation, self._scale, self._translation)[0] + fieldassignment = self._fitter._modelCoordinatesField.createFieldassignment(modelCoordinatesTransformed) + result = fieldassignment.assign() + assert result in [ RESULT_OK, RESULT_WARNING_PART_DONE ], "Align: Failed to transform model" + self._fitter.updateModelReferenceCoordinates() + del fieldassignment + del modelCoordinatesTransformed + self._fitter.calculateDataProjections() + self.setHasRun(True) + + def _doAlignMarkers(self): + """ + Prepare and invoke alignment to markers. + """ + fieldmodule = self._fitter._fieldmodule + markerGroup = self._fitter.getMarkerGroup() + assert markerGroup, "Align: No marker group to align with" + markerPrefix = markerGroup.getName() + modelCoordinates = self._fitter.getModelCoordinatesField() + componentsCount = modelCoordinates.getNumberOfComponents() + + markerNodeGroup, markerLocation, markerCoordinates, markerName = self._fitter.getMarkerModelFields() + assert markerNodeGroup and markerCoordinates and markerName, "Align: No marker group, coordinates or name fields" + modelMarkers = getNodeNameCentres(markerNodeGroup, markerCoordinates, markerName) + + markerDataGroup, markerDataCoordinates, markerDataName = self._fitter.getMarkerDataFields() + assert markerDataGroup and markerDataCoordinates and markerDataName, "Align: No marker data group, coordinates or name fields" + dataMarkers = getNodeNameCentres(markerDataGroup, markerDataCoordinates, markerDataName) + + # match model and data markers, warn of missing markers + markerMap = {} + for name, modelx in modelMarkers.items(): + datax = dataMarkers.get(name) + if datax: + markerMap[name] = ( modelx, datax ) + if self.getDiagnosticLevel() > 0: + for name in modelMarkers: + datax = dataMarkers.get(name) + if datax: + print("Align: Found marker " + name + " in model and data") + for name in modelMarkers: + if not markerMap.get(name): + print("Align: Model marker " + name + " not found in data") + for name in dataMarkers: + if not markerMap.get(name): + print("Align: Data marker " + name + " not found in model") + + self._optimiseAlignment(markerMap) + + def _optimiseAlignment(self, markerMap): + """ + Calculate transformation from modelCoordinates to dataMarkers + over the markers, by scaling, translating and rotating model. + On success, sets transformation parameters in object. + :param markerMap: dict name -> (modelCoordinates, dataCoordinates) + """ + assert len(markerMap) >= 3, "Align: Only " + str(len(markerMap)) + " markers - need at least 3" + region = self._fitter._context.createRegion() + fieldmodule = region.getFieldmodule() + with ChangeManager(fieldmodule): + modelCoordinates = fieldmodule.createFieldFiniteElement(3) + dataCoordinates = fieldmodule.createFieldFiniteElement(3) + nodes = fieldmodule.findNodesetByFieldDomainType(Field.DOMAIN_TYPE_NODES) + nodetemplate = nodes.createNodetemplate() + nodetemplate.defineField(modelCoordinates) + nodetemplate.defineField(dataCoordinates) + fieldcache = fieldmodule.createFieldcache() + for name, positions in markerMap.items(): + modelx = positions[0] + datax = positions[1] + node = nodes.createNode(-1, nodetemplate) + fieldcache.setNode(node) + result1 = modelCoordinates.assignReal(fieldcache, positions[0]) + result2 = dataCoordinates.assignReal(fieldcache, positions[1]) + assert (result1 == RESULT_OK) and (result2 == RESULT_OK), "Align: Failed to set up data for alignment to markers optimisation" + del fieldcache + modelCoordinatesTransformed, rotation, scale, translation = createFieldsTransformations(modelCoordinates) + # create objective = sum of squares of vector from modelCoordinatesTransformed to dataCoordinates + markerDiff = fieldmodule.createFieldSubtract(dataCoordinates, modelCoordinatesTransformed) + objective = fieldmodule.createFieldNodesetSumSquares(markerDiff, nodes) + assert objective.isValid(), "Align: Failed to set up objective function for alignment to markers optimisation" + + # future: pre-fit to avoid gimbal lock + + optimisation = fieldmodule.createOptimisation() + optimisation.setMethod(Optimisation.METHOD_LEAST_SQUARES_QUASI_NEWTON) + optimisation.addObjectiveField(objective) + optimisation.addIndependentField(rotation) + optimisation.addIndependentField(scale) + optimisation.addIndependentField(translation) + + result = optimisation.optimise() + if self.getDiagnosticLevel() > 1: + solutionReport = optimisation.getSolutionReport() + print(solutionReport) + assert result == RESULT_OK, "Align: Alignment to markers optimisation failed" + + fieldcache = fieldmodule.createFieldcache() + result1, self._rotation = rotation.evaluateReal(fieldcache, 3) + result2, self._scale = scale.evaluateReal(fieldcache, 1) + result3, self._translation = translation.evaluateReal(fieldcache, 3) + assert (result1 == RESULT_OK) and (result2 == RESULT_OK) and (result3 == RESULT_OK), "Align: Failed to evaluate transformation for alignment to markers" diff --git a/src/scaffoldfitter/fitterstepfit.py b/src/scaffoldfitter/fitterstepfit.py new file mode 100644 index 0000000..57616c4 --- /dev/null +++ b/src/scaffoldfitter/fitterstepfit.py @@ -0,0 +1,237 @@ +""" +Fit step for gross alignment and scale. +""" + +from opencmiss.utils.zinc.field import assignFieldParameters, createFieldsDisplacementGradients +from opencmiss.utils.zinc.general import ChangeManager +from opencmiss.zinc.field import Field, FieldFindMeshLocation +from opencmiss.zinc.optimisation import Optimisation +from opencmiss.zinc.result import RESULT_OK +from scaffoldfitter.fitter import Fitter, FitterStep + +class FitterStepFit(FitterStep): + + _jsonTypeId = "_FitterStepFit" + + def __init__(self, fitter : Fitter): + super(FitterStepFit, self).__init__(fitter) + self._markerWeight = 1.0 + self._strainPenaltyWeight = 0.0 + self._curvaturePenaltyWeight = 0.0 + self._edgeDiscontinuityPenaltyWeight = 0.0 + self._numberOfIterations = 1 + self._updateReferenceState = False + + @classmethod + def getJsonTypeId(cls): + return cls._jsonTypeId + + def decodeSettingsJSONDict(self, dct : dict): + """ + Decode definition of step from JSON dict. + """ + assert self._jsonTypeId in dct + self._markerWeight = dct["markerWeight"] + self._strainPenaltyWeight = dct["strainPenaltyWeight"] + self._curvaturePenaltyWeight = dct["curvaturePenaltyWeight"] + self._edgeDiscontinuityPenaltyWeight = dct["edgeDiscontinuityPenaltyWeight"] + self._numberOfIterations = dct["numberOfIterations"] + self._updateReferenceState = dct["updateReferenceState"] + + def encodeSettingsJSONDict(self) -> dict: + """ + Encode definition of step in dict. + :return: Settings in a dict ready for passing to json.dump. + """ + return { + self._jsonTypeId : True, + "markerWeight" : self._markerWeight, + "strainPenaltyWeight" : self._strainPenaltyWeight, + "curvaturePenaltyWeight" : self._curvaturePenaltyWeight, + "edgeDiscontinuityPenaltyWeight" : self._edgeDiscontinuityPenaltyWeight, + "numberOfIterations" : self._numberOfIterations, + "updateReferenceState" : self._updateReferenceState + } + + def getMarkerWeight(self): + return self._markerWeight + + def setMarkerWeight(self, weight): + assert weight >= 0.0 + self._markerWeight = weight + + def getStrainPenaltyWeight(self): + return self._strainPenaltyWeight + + def setStrainPenaltyWeight(self, weight): + assert weight >= 0.0 + self._strainPenaltyWeight = weight + + def getCurvaturePenaltyWeight(self): + return self._curvaturePenaltyWeight + + def setCurvaturePenaltyWeight(self, weight): + assert weight >= 0.0 + self._curvaturePenaltyWeight = weight + + def getEdgeDiscontinuityPenaltyWeight(self): + return self._edgeDiscontinuityPenaltyWeight + + def setEdgeDiscontinuityPenaltyWeight(self, weight): + assert weight >= 0.0 + self._edgeDiscontinuityPenaltyWeight = weight + + def getNumberOfIterations(self): + return self._numberOfIterations + + def setNumberOfIterations(self, numberOfIterations): + assert numberOfIterations >= 0 # 0 iterations performs projections only + self._numberOfIterations = numberOfIterations + + def isUpdateReferenceState(self): + return self._updateReferenceState + + def setUpdateReferenceState(self, updateReferenceState): + self._updateReferenceState = updateReferenceState + + def run(self): + """ + Fit model geometry parameters to data. + """ + fieldmodule = self._fitter._region.getFieldmodule() + optimisation = fieldmodule.createOptimisation() + optimisation.setMethod(Optimisation.METHOD_LEAST_SQUARES_QUASI_NEWTON) + optimisation.addIndependentField(self._fitter.getModelCoordinatesField()) + optimisation.setAttributeInteger(Optimisation.ATTRIBUTE_MAXIMUM_ITERATIONS, 5) + + dataProjectionObjective = [ None, None ] + dataProjectionObjectiveComponentsCount = [ 0, 0 ] + markerObjectiveField = None + deformationPenaltyObjective = None + edgeDiscontinuityPenaltyObjective = None + for dimension in range(1, 3): + if self._fitter.getDataProjectionNodesetGroup(dimension).getSize() > 0: + dataProjectionObjective[dimension - 1] = self.createDataProjectionObjectiveField(dimension) + dataProjectionObjectiveComponentsCount[dimension - 1] = dataProjectionObjective[dimension - 1].getNumberOfComponents() + result = optimisation.addObjectiveField(dataProjectionObjective[dimension - 1]) + assert result == RESULT_OK, "Fit Geometry: Could not add data projection objective field for dimension " + str(dimension) + if self._fitter.getMarkerGroup() and (self._fitter.getMarkerDataLocationNodesetGroup().getSize() > 0) and (self._markerWeight > 0.0): + markerObjectiveField = self.createMarkerObjectiveField(self._markerWeight) + result = optimisation.addObjectiveField(markerObjectiveField) + assert result == RESULT_OK, "Fit Geometry: Could not add marker objective field" + if (self._strainPenaltyWeight > 0.0) or (self._curvaturePenaltyWeight > 0.0): + deformationPenaltyObjective = self.createDeformationPenaltyObjectiveField() + result = optimisation.addObjectiveField(deformationPenaltyObjective) + assert result == RESULT_OK, "Fit Geometry: Could not add strain/curvature penalty objective field" + if self._edgeDiscontinuityPenaltyWeight > 0.0: + edgeDiscontinuityPenaltyObjective = self.createEdgeDiscontinuityPenaltyObjectiveField() + result = optimisation.addObjectiveField(edgeDiscontinuityPenaltyObjective) + assert result == RESULT_OK, "Fit Geometry: Could not add edge discontinuity penalty objective field" + + fieldcache = fieldmodule.createFieldcache() + for iter in range(self._numberOfIterations): + if self.getDiagnosticLevel() > 0: + print("-------- Iteration", iter + 1) + for d in range(2): + if dataProjectionObjective[d]: + result, objective = dataProjectionObjective[d].evaluateReal(fieldcache, dataProjectionObjectiveComponentsCount[d]) + if self.getDiagnosticLevel() > 0: + print(" " + str(d + 1) + "-D data projection objective", objective) + result = optimisation.optimise() + if self.getDiagnosticLevel() > 1: + solutionReport = optimisation.getSolutionReport() + print(solutionReport) + assert result == RESULT_OK, "Fit Geometry: Optimisation failed with result " + str(result) + self._fitter.calculateDataProjections() + if self.getDiagnosticLevel() > 0: + print("--------") + + for d in range(2): + if dataProjectionObjective[d]: + result, objective = dataProjectionObjective[d].evaluateReal(fieldcache, dataProjectionObjectiveComponentsCount[d]) + if self.getDiagnosticLevel() > 0: + print("END " + str(d + 1) + "-D data projection objective", objective) + + if self._updateReferenceState: + self._fitter.updateModelReferenceCoordinates() + + self.setHasRun(True) + + def createDataProjectionObjectiveField(self, dimension): + """ + Get FieldNodesetSumSquares objective for data projected onto mesh of dimension. + Minimises length in projection direction, allowing sliding fit. + Only call if self._fitter.getDataProjectionNodesetGroup().getSize() > 0 + :param dimension: Mesh dimension 1 or 2. + :return: Zinc FieldNodesetSumSquares. + """ + fieldmodule = self._fitter.getFieldmodule() + with ChangeManager(fieldmodule): + dataProjectionDelta = self._fitter.getDataProjectionDeltaField(dimension) + #dataProjectionInDirection = fieldmodule.createFieldDotProduct(dataProjectionDelta, self._fitter.getDataProjectionDirectionField()) + #dataProjectionInDirection = fieldmodule.createFieldMagnitude(dataProjectionDelta) + dataProjectionInDirection = dataProjectionDelta + dataProjectionObjective = fieldmodule.createFieldNodesetSumSquares(dataProjectionInDirection, self._fitter.getDataProjectionNodesetGroup(dimension)) + return dataProjectionObjective + + def createMarkerObjectiveField(self, weight): + """ + Only call if self._fitter.getMarkerGroup() and (self._fitter.getMarkerDataLocationNodesetGroup().getSize() > 0) and (self._markerWeight > 0.0) + For marker datapoints with locations in model, creates a FieldNodesetSumSquares + of coordinate projections to those locations. + :return: Zinc FieldNodesetSumSquares. + """ + fieldmodule = self._fitter.getFieldmodule() + with ChangeManager(fieldmodule): + markerDataLocation, markerDataLocationCoordinates, markerDataDelta = self._fitter.getMarkerDataLocationFields() + markerDataWeightedDelta = markerDataDelta*fieldmodule.createFieldConstant([ weight ]*markerDataDelta.getNumberOfComponents()) + markerDataObjective = fieldmodule.createFieldNodesetSumSquares(markerDataWeightedDelta, self._fitter.getMarkerDataLocationNodesetGroup()) + return markerDataObjective + + def createDeformationPenaltyObjectiveField(self): + """ + Only call if (self._strainPenaltyWeight > 0.0) or (self._curvaturePenaltyWeight > 0.0) + :return: Zinc FieldMeshIntegralSquares, or None if not weighted. + """ + numberOfGaussPoints = 3 + fieldmodule = self._fitter.getFieldmodule() + mesh = self._fitter.getHighestDimensionMesh() + with ChangeManager(fieldmodule): + displacementGradient1, displacementGradient2 = createFieldsDisplacementGradients(self._fitter.getModelCoordinatesField(), self._fitter.getModelReferenceCoordinatesField(), mesh) + if self._strainPenaltyWeight > 0.0: + weightedDisplacementGradient1 = displacementGradient1*fieldmodule.createFieldConstant([ self._strainPenaltyWeight ]*displacementGradient1.getNumberOfComponents()) + else: + weightedDisplacementGradient1 = None + if self._curvaturePenaltyWeight > 0.0: + weightedDisplacementGradient2 = displacementGradient2*fieldmodule.createFieldConstant([ self._curvaturePenaltyWeight ]*displacementGradient2.getNumberOfComponents()) + else: + weightedDisplacementGradient2 = None + + if weightedDisplacementGradient1: + if weightedDisplacementGradient2: + deformationField = fieldmodule.createFieldConcatenate([ weightedDisplacementGradient1, weightedDisplacementGradient2 ]) + else: + deformationField = weightedDisplacementGradient1 + elif weightedDisplacementGradient2: + deformationField = weightedDisplacementGradient2 + else: + return None + + deformationPenaltyObjective = fieldmodule.createFieldMeshIntegralSquares(deformationField, self._fitter.getModelReferenceCoordinatesField(), mesh) + deformationPenaltyObjective.setNumbersOfPoints(numberOfGaussPoints) + return deformationPenaltyObjective + + def createEdgeDiscontinuityPenaltyObjectiveField(self): + """ + Only call if self._edgeDiscontinuityPenaltyWeight > 0.0 + :return: Zinc FieldMeshIntegralSquares, or None if not weighted. + """ + numberOfGaussPoints = 3 + fieldmodule = self._fitter.getFieldmodule() + lineMesh = fieldmodule.findMeshByDimension(1) + with ChangeManager(fieldmodule): + edgeDiscontinuity = fieldmodule.createFieldEdgeDiscontinuity(self._fitter.getModelCoordinatesField()) + weightedEdgeDiscontinuity = edgeDiscontinuity*fieldmodule.createFieldConstant(self._edgeDiscontinuityPenaltyWeight) + edgeDiscontinuityPenaltyObjective = fieldmodule.createFieldMeshIntegralSquares(weightedEdgeDiscontinuity, self._fitter.getModelReferenceCoordinatesField(), lineMesh) + edgeDiscontinuityPenaltyObjective.setNumbersOfPoints(numberOfGaussPoints) + return edgeDiscontinuityPenaltyObjective diff --git a/tests/resources/cube_to_sphere.exf b/tests/resources/cube_to_sphere.exf new file mode 100644 index 0000000..cd78694 --- /dev/null +++ b/tests/resources/cube_to_sphere.exf @@ -0,0 +1,400 @@ +EX Version: 2 +Region: / +!#nodeset nodes +Shape. Dimension=0 +#Fields=1 +1) coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=4 (value,d/ds1,d/ds2,d/ds3) + y. #Values=4 (value,d/ds1,d/ds2,d/ds3) + z. #Values=4 (value,d/ds1,d/ds2,d/ds3) +Node: 1 + 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 2 + 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 3 + 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 4 + 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 5 + 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 6 + 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 7 + 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Node: 8 + 1.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 + 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +Shape. Dimension=0 +#Fields=2 +1) marker_location, field, rectangular cartesian, element_xi, #Components=1, host mesh=mesh3d, host mesh dimension=3 + 1. #Values=1 (value) +2) marker_name, field, rectangular cartesian, string, #Components=1 + 1. #Values=1 (value) +Node: 1000001 + 1 5.000000000000000e-01 5.000000000000000e-01 0.000000000000000e+00 + bottom +Node: 1000002 + 1 0.000000000000000e+00 0.000000000000000e+00 5.000000000000000e-01 + left +Node: 1000003 + 1 1.000000000000000e+00 1.000000000000000e+00 5.000000000000000e-01 + right +Node: 1000004 + 1 5.000000000000000e-01 5.000000000000000e-01 1.000000000000000e+00 + top +!#mesh mesh1d, dimension=1, nodeset=nodes +Shape. Dimension=1, line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 1 +Element: 2 +Element: 3 +Element: 4 +Element: 5 +Element: 6 +Element: 7 +Element: 8 +Element: 9 +Element: 10 +Element: 11 +Element: 12 +!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes +Shape. Dimension=2, line*line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 1 + Faces: + 1 2 3 4 +Element: 2 + Faces: + 5 6 7 8 +Element: 3 + Faces: + 9 10 1 5 +Element: 4 + Faces: + 11 12 2 6 +Element: 5 + Faces: + 3 7 9 11 +Element: 6 + Faces: + 4 8 10 12 +!#mesh mesh3d, dimension=3, face mesh=mesh2d, nodeset=nodes +Shape. Dimension=3, line*line*line +#Scale factor sets=0 +#Nodes=8 +#Fields=1 +1) coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. c.Hermite*c.Hermite*c.Hermite, no modify, standard node based. + #Nodes=8 + 1. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 1. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 2. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 2. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 3. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 3. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 4. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 4. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 5. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 5. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 6. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 6. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 7. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 7. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 8. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 8. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + y. c.Hermite*c.Hermite*c.Hermite, no modify, standard node based. + #Nodes=8 + 1. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 1. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 2. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 2. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 3. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 3. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 4. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 4. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 5. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 5. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 6. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 6. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 7. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 7. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 8. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 8. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + z. c.Hermite*c.Hermite*c.Hermite, no modify, standard node based. + #Nodes=8 + 1. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 1. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 2. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 2. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 3. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 3. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 4. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 4. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 5. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 5. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 6. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 6. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 7. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 7. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero + 8. #Values=3 + Value labels: value d/ds1 d/ds2 + 0. #Values=1 + Value labels: zero + 8. #Values=1 + Value labels: d/ds3 + 0. #Values=3 + Value labels: zero zero zero +Element: 1 + Faces: + 1 2 3 4 5 6 + Nodes: + 1 2 3 4 5 6 7 8 + Group name: bottom +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 1 +Node: 2 +Node: 3 +Node: 4 +!#mesh mesh1d, dimension=1, nodeset=nodes +Shape. Dimension=1, line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 3 +Element: 7 +Element: 9 +Element: 11 +!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes +Shape. Dimension=2, line*line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 5 + Group name: marker +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 1000001 +Node: 1000002 +Node: 1000003 +Node: 1000004 + Group name: sides +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 1 +Node: 2 +Node: 3 +Node: 4 +Node: 5 +Node: 6 +Node: 7 +Node: 8 +!#mesh mesh1d, dimension=1, nodeset=nodes +Shape. Dimension=1, line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 1 +Element: 2 +Element: 3 +Element: 4 +Element: 5 +Element: 6 +Element: 7 +Element: 8 +Element: 9 +Element: 10 +Element: 11 +Element: 12 +!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes +Shape. Dimension=2, line*line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 1 +Element: 2 +Element: 3 +Element: 4 + Group name: top +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 5 +Node: 6 +Node: 7 +Node: 8 +!#mesh mesh1d, dimension=1, nodeset=nodes +Shape. Dimension=1, line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 4 +Element: 8 +Element: 10 +Element: 12 +!#mesh mesh2d, dimension=2, face mesh=mesh1d, nodeset=nodes +Shape. Dimension=2, line*line +#Scale factor sets=0 +#Nodes=0 +#Fields=0 +Element: 6 diff --git a/tests/resources/cube_to_sphere_data_random.exf b/tests/resources/cube_to_sphere_data_random.exf new file mode 100644 index 0000000..df35738 --- /dev/null +++ b/tests/resources/cube_to_sphere_data_random.exf @@ -0,0 +1,866 @@ +EX Version: 2 +Region: / +!#nodeset datapoints +Shape. Dimension=0 +#Fields=1 +1) data_coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=1 (value) + y. #Values=1 (value) + z. #Values=1 (value) +Node: 1 + 8.127885782815432e-02 + 1.384642158246737e-01 + -4.730473585490702e-01 +Node: 2 + -6.124387617180985e-02 + 5.279891597195788e-02 + -4.931691876344948e-01 +Node: 3 + -1.709978132375876e-01 + 6.990974695120276e-02 + -4.642635806774232e-01 +Node: 4 + -1.827908530010493e-01 + 1.155753355184168e-02 + -4.647528303441892e-01 +Node: 5 + 3.442951169644994e-02 + -1.564746590768734e-01 + -4.734186660038504e-01 +Node: 6 + 6.295815555042669e-02 + -1.659330886146003e-01 + -4.673392997663131e-01 +Node: 7 + 1.274553989427585e-02 + -2.446596387837532e-02 + -4.987761715693815e-01 +Node: 8 + 1.635281205205502e-01 + -6.303297746524651e-02 + -4.682471805082765e-01 +Node: 9 + 1.752616315425237e-01 + -1.331047692695432e-02 + -4.678509035196187e-01 +Node: 10 + 3.311868119702572e-01 + 1.146261360257874e-01 + -3.556745197079050e-01 +Node: 11 + 1.278737743381185e-01 + 3.191410653918657e-01 + -3.617388813187926e-01 +Node: 12 + 1.347968667494562e-01 + 3.993643493285188e-01 + -2.676752569343143e-01 +Node: 13 + 7.147713652042686e-02 + 2.750609094765267e-01 + -4.106300672316757e-01 +Node: 14 + -3.890260803579240e-03 + 2.936256618461847e-01 + -4.041233251933907e-01 +Node: 15 + -2.781866954901989e-03 + 3.576762586414576e-01 + -3.487698646719835e-01 +Node: 16 + -3.764801439650819e-01 + 1.077824277316650e-01 + -3.099156735110160e-01 +Node: 17 + -3.297655302613394e-01 + 4.026011828730260e-02 + -3.723992012543707e-01 +Node: 18 + -3.367606897178697e-01 + -1.793664064232528e-02 + -3.679903923681171e-01 +Node: 19 + -3.074453711004010e-01 + -6.051675979995395e-02 + -3.883683450569380e-01 +Node: 20 + -4.063429183305032e-01 + -1.328356705607575e-01 + -2.580151320646145e-01 +Node: 21 + -2.729131834447396e-01 + -1.048009171013258e-01 + -4.046856263333486e-01 +Node: 22 + -2.973794052159342e-01 + -1.565910095918145e-01 + -3.691024919157517e-01 +Node: 23 + -3.338882772627539e-01 + -1.562405854525425e-01 + -3.368792627800896e-01 +Node: 24 + -9.920809414029344e-02 + -3.123464052470948e-01 + -3.761628844321663e-01 +Node: 25 + -4.216388144206367e-02 + -2.925314324148138e-01 + -4.021930013666090e-01 +Node: 26 + 2.023970451140369e-02 + -2.657119209644258e-01 + -4.223726416376690e-01 +Node: 27 + 1.899338151349948e-02 + -4.124612531750902e-01 + -2.809586299401809e-01 +Node: 28 + 3.554662190778812e-02 + -2.708920352851031e-01 + -4.178608928322446e-01 +Node: 29 + 1.240387584271821e-01 + -3.467025110765755e-01 + -3.371215677604788e-01 +Node: 30 + 1.995112426483726e-01 + -2.513704904995668e-01 + -3.825997454060966e-01 +Node: 31 + 2.128170646439904e-01 + -3.002947878003736e-01 + -3.373972264628650e-01 +Node: 32 + 2.387922119787972e-01 + -1.696997889039514e-01 + -4.044689821712003e-01 +Node: 33 + 3.395145877722356e-01 + -2.676188272932835e-01 + -2.504376174583788e-01 +Node: 34 + 3.106729836686143e-01 + -2.901636736928639e-03 + -3.910548949478478e-01 +Node: 35 + 4.980235866183390e-01 + 8.971728229307512e-03 + -3.548505091110003e-02 +Node: 36 + 4.673984874477937e-01 + 1.093235568383570e-01 + -1.351246625993578e-01 +Node: 37 + 4.669866545572816e-01 + 1.243043044047040e-01 + -1.247245075514978e-01 +Node: 38 + 4.318430097600009e-01 + 1.444855501243334e-01 + -2.036186844218722e-01 +Node: 39 + 4.372151638567493e-01 + 1.550070048142539e-01 + -1.834340308657512e-01 +Node: 40 + 4.673601056364388e-01 + 1.741171408896335e-01 + -2.256543933003704e-02 +Node: 41 + 2.451931401283649e-01 + 4.294558919659939e-01 + -6.949826812686932e-02 +Node: 42 + 1.916480372024419e-01 + 4.580248116139006e-01 + -5.080432589800116e-02 +Node: 43 + 1.884368259622198e-01 + 4.607919329251485e-01 + -3.586003346213541e-02 +Node: 44 + 3.167887988599520e-02 + 4.461370122635714e-01 + -2.218684309356638e-01 +Node: 45 + -1.898567180858996e-02 + 4.480776067453792e-01 + -2.193232555994346e-01 +Node: 46 + -8.130980260374818e-02 + 4.805619852190419e-01 + -1.066882429539814e-01 +Node: 47 + -2.562217969051391e-01 + 3.685743839088007e-01 + -2.191023981572653e-01 +Node: 48 + -3.117412880518273e-01 + 3.610571411430015e-01 + -1.453039181573705e-01 +Node: 49 + -4.278137578927707e-01 + 2.089434360598872e-01 + -1.494482477247091e-01 +Node: 50 + -4.848167675611508e-01 + 1.196412637169922e-01 + -4.574255232488424e-03 +Node: 51 + -4.846561449927872e-01 + 2.279513069598848e-02 + -1.165341202490377e-01 +Node: 52 + -4.724446714441017e-01 + -3.644245634989011e-02 + -1.569242932187317e-01 +Node: 53 + -4.687422941555585e-01 + -1.438918320806760e-01 + -9.188369233730231e-02 +Node: 54 + -4.212907548799074e-01 + -2.678728639510137e-01 + -5.986423562983369e-03 +Node: 55 + -4.281521659292933e-01 + -2.572315096933382e-01 + -4.634092873611090e-03 +Node: 56 + -3.590380027474099e-01 + -3.103916812856942e-01 + -1.528853494829155e-01 +Node: 57 + -3.092468769485129e-01 + -3.629950970346725e-01 + -1.461495442271767e-01 +Node: 58 + -2.916704814573099e-01 + -3.489990739323678e-01 + -2.059586131756811e-01 +Node: 59 + -2.740009877638948e-01 + -3.762129072505684e-01 + -1.800262481906458e-01 +Node: 60 + -2.115415924376781e-01 + -3.885501915301395e-01 + -2.312519195882584e-01 +Node: 61 + -2.200559694534270e-01 + -4.472266303824269e-01 + -2.665966370818576e-02 +Node: 62 + -1.968677996680498e-01 + -4.235546305917983e-01 + -1.768091537157719e-01 +Node: 63 + -6.697712735941926e-02 + -4.566586681841989e-01 + -1.889172393072476e-01 +Node: 64 + -5.030080476370402e-02 + -4.740387037513304e-01 + -1.470434192318968e-01 +Node: 65 + 1.231844585022696e-02 + -4.937020008685217e-01 + -7.217278964704971e-02 +Node: 66 + 4.183122107257346e-02 + -4.699460984105852e-01 + -1.635585217693885e-01 +Node: 67 + 9.147092389834681e-02 + -4.288456402308599e-01 + -2.382839260624621e-01 +Node: 68 + 1.039736532280831e-01 + -4.863557605215594e-01 + -4.418597368032735e-02 +Node: 69 + 1.557267201375375e-01 + -4.657964514329755e-01 + -8.866367334071400e-02 +Node: 70 + 1.631424197855549e-01 + -4.604299897371506e-01 + -1.014852744058217e-01 +Node: 71 + 1.986652953829841e-01 + -4.290280841739118e-01 + -1.603499675976163e-01 +Node: 72 + 2.374768138975956e-01 + -3.922389998232095e-01 + -1.973616509613717e-01 +Node: 73 + 2.401226515855199e-01 + -4.039517471314233e-01 + -1.696848667946363e-01 +Node: 74 + 2.889732902591439e-01 + -3.475585615234364e-01 + -2.124825801713883e-01 +Node: 75 + 3.596111302442443e-01 + -3.446160054484798e-01 + -3.180304139401858e-02 +Node: 76 + 3.642093705680586e-01 + -3.416643753117098e-01 + -1.854966874802768e-03 +Node: 77 + 3.924777074403626e-01 + -3.083234263662635e-01 + -8.765549561791694e-03 +Node: 78 + 4.301822463812039e-01 + -1.310042275205833e-01 + -2.161733796728808e-01 +Node: 79 + 4.973596948203152e-01 + -3.080203718001219e-02 + -2.853856563943751e-02 +Node: 80 + 4.960752655358777e-01 + 2.937161932596009e-02 + 4.864438917943417e-02 +Node: 81 + 4.542146012544901e-01 + 1.834402596011540e-03 + 2.084474017503158e-01 +Node: 82 + 4.907003725714129e-01 + 5.358334955143060e-02 + 7.336557082822537e-02 +Node: 83 + 3.629100932458149e-01 + 3.417801723103828e-01 + 2.081020185710540e-02 +Node: 84 + 3.109406174005950e-01 + 3.527121996862843e-01 + 1.664064786381203e-01 +Node: 85 + 2.152836087759273e-01 + 3.780399514388081e-01 + 2.458745719492117e-01 +Node: 86 + 1.890435669656770e-01 + 4.230615263379308e-01 + 1.850788703735548e-01 +Node: 87 + 1.710933935652429e-01 + 4.232175297846191e-01 + 2.011975594021583e-01 +Node: 88 + 2.028385867411695e-02 + 4.379978340398688e-01 + 2.386772176573412e-01 +Node: 89 + -4.023477579430475e-02 + 4.355909077153162e-01 + 2.411432622522584e-01 +Node: 90 + -2.161059726337622e-01 + 3.961968055789579e-01 + 2.140781866118735e-01 +Node: 91 + -2.570753266396119e-01 + 4.227094927618333e-01 + 6.576950748748610e-02 +Node: 92 + -2.939703214052857e-01 + 3.666295265889122e-01 + 1.678641043560445e-01 +Node: 93 + -3.017025520846511e-01 + 3.402645332609093e-01 + 2.050999011514126e-01 +Node: 94 + -2.912412366015534e-01 + 3.228785904872973e-01 + 2.454848235218541e-01 +Node: 95 + -4.049716328482794e-01 + 2.924009161921814e-01 + 4.310969611548585e-03 +Node: 96 + -3.940577898004557e-01 + 2.517075924945823e-01 + 1.749699352683774e-01 +Node: 97 + -4.072133929486001e-01 + 2.063624179820377e-01 + 2.018446690364924e-01 +Node: 98 + -4.696287302061073e-01 + 1.375597101536553e-01 + 9.651360665725556e-02 +Node: 99 + -4.946817387351657e-01 + 2.806511503899326e-02 + 5.943087829017649e-02 +Node: 100 + -4.939539154093525e-01 + 2.198834817275711e-02 + 6.690687184387263e-02 +Node: 101 + -4.856012800429474e-01 + -8.173462284398166e-03 + 1.153854129072362e-01 +Node: 102 + -4.790933453524106e-01 + -6.765388192896370e-03 + 1.403545932022585e-01 +Node: 103 + -4.573771964706223e-01 + -8.342456437605733e-02 + 1.809981428886001e-01 +Node: 104 + -4.551487080674403e-01 + -1.355398052982584e-01 + 1.513286864775364e-01 +Node: 105 + -3.869608305975607e-01 + -2.919513364806946e-01 + 1.176308498830654e-01 +Node: 106 + -2.509506369452381e-01 + -3.707960569928038e-01 + 2.210027170045938e-01 +Node: 107 + -1.195879233192418e-01 + -4.798344626195533e-01 + 6.576027135002124e-02 +Node: 108 + -2.647694083790244e-02 + -4.654481716597961e-01 + 1.783727969576704e-01 +Node: 109 + 8.945459117034853e-02 + -4.639581573323459e-01 + 1.594827564779802e-01 +Node: 110 + 1.731602301899047e-01 + -4.608614036175178e-01 + 8.320337667650607e-02 +Node: 111 + 2.299075181437582e-01 + -3.937572128921825e-01 + 2.041809070509419e-01 +Node: 112 + 3.331683252216983e-01 + -3.576232825972670e-01 + 9.878893081505465e-02 +Node: 113 + 3.944812519087850e-01 + -2.395315355034835e-01 + 1.900958154093238e-01 +Node: 114 + 3.909974119624821e-01 + -2.207311515339047e-01 + 2.188803596592555e-01 +Node: 115 + 4.652295167172116e-01 + -1.790460872027851e-01 + 2.457598407175861e-02 +Node: 116 + 4.065530303469233e-01 + -1.556031534727774e-01 + 2.448742176523272e-01 +Node: 117 + 4.028092463117576e-01 + 1.239647130696859e-01 + 2.676393498608062e-01 +Node: 118 + 3.488441421423396e-01 + 9.643343221618472e-02 + 3.437981981107989e-01 +Node: 119 + 3.455694986471313e-01 + 1.483320353365281e-01 + 3.284368471475808e-01 +Node: 120 + 2.691224402557247e-01 + 1.127863793095867e-01 + 4.051108006440815e-01 +Node: 121 + 2.850538463534611e-01 + 1.620311468892398e-01 + 3.768646948263047e-01 +Node: 122 + 1.994228258373812e-01 + 1.704860626837192e-01 + 4.248826688213912e-01 +Node: 123 + 2.944606964360130e-01 + 3.088818831794425e-01 + 2.593201297674209e-01 +Node: 124 + 2.528871480059128e-01 + 2.889837094928491e-01 + 3.196711759171509e-01 +Node: 125 + 1.564936361674638e-01 + 2.187163608730391e-01 + 4.207532230724722e-01 +Node: 126 + 2.383000906003741e-01 + 3.597381714544207e-01 + 2.518009641556357e-01 +Node: 127 + 1.735646593916146e-01 + 2.949684163083177e-01 + 3.636706594165042e-01 +Node: 128 + 1.548774214537903e-01 + 2.362021621314739e-01 + 4.120748005051359e-01 +Node: 129 + 5.966395008347802e-02 + 2.536259562934027e-01 + 4.259808952046855e-01 +Node: 130 + 1.376327192294579e-02 + 4.166933360522393e-01 + 2.746704036920920e-01 +Node: 131 + -7.385987246522578e-02 + 3.460756690342382e-01 + 3.522621955634233e-01 +Node: 132 + -7.554251576512647e-02 + 2.622193655832885e-01 + 4.179618380213219e-01 +Node: 133 + -1.625615003720760e-01 + 3.870379475451885e-01 + 2.702313917296651e-01 +Node: 134 + -2.418989646224502e-01 + 3.339436220590489e-01 + 2.821435633659450e-01 +Node: 135 + -2.767166350672730e-01 + 1.174912030667964e-01 + 3.983404677266091e-01 +Node: 136 + -2.722696933950944e-01 + 6.225098735514816e-02 + 4.137436177798811e-01 +Node: 137 + -3.468031292177481e-01 + -5.093351147881771e-02 + 3.556484267573538e-01 +Node: 138 + -3.078462603101871e-01 + -1.357557813372600e-01 + 3.685815326616413e-01 +Node: 139 + -3.073314257293769e-01 + -1.795252339436544e-01 + 3.507436202207410e-01 +Node: 140 + -3.141435304530485e-01 + -2.595034133488634e-01 + 2.895405375641339e-01 +Node: 141 + -2.884768693452279e-01 + -2.265516643405277e-01 + 3.384802414359407e-01 +Node: 142 + -1.737652597012954e-01 + -1.921049002287978e-01 + 4.269574146224257e-01 +Node: 143 + -1.925862084666152e-01 + -3.438437939586975e-01 + 3.066205800636486e-01 +Node: 144 + -9.913027997753365e-02 + -3.989316487840204e-01 + 2.843552146748285e-01 +Node: 145 + -9.680737434037784e-02 + -3.902373686944445e-01 + 2.964050518031902e-01 +Node: 146 + -9.845683778827925e-02 + -3.680424027877757e-01 + 3.234740095128072e-01 +Node: 147 + 5.161349495347718e-03 + -3.702409864745316e-01 + 3.349692716998555e-01 +Node: 148 + 3.370693715890486e-01 + -2.071807864123087e-01 + 3.044227959935665e-01 +Node: 149 + 3.820549802484084e-01 + -1.721800108798260e-01 + 2.716183024439118e-01 +Node: 150 + 2.822145622097533e-01 + -6.043880892732824e-02 + 4.073380973672912e-01 +Node: 151 + 2.151309409541949e-01 + 3.111991867854853e-02 + 4.500776213607793e-01 +Node: 152 + 1.427907685525803e-01 + 5.843053723866751e-02 + 4.751616927103051e-01 +Node: 153 + 1.287567740004467e-01 + 1.965596597803942e-01 + 4.407527391711622e-01 +Node: 154 + 1.459260468371241e-02 + 1.236794055987093e-01 + 4.840112735586134e-01 +Node: 155 + 1.367143842438823e-02 + 1.909521620067020e-01 + 4.613595816169517e-01 +Node: 156 + -2.551009261342680e-02 + 8.317726984755160e-02 + 4.923603851912634e-01 +Node: 157 + -4.350107679589271e-02 + 4.553643823738444e-02 + 4.955297418419164e-01 +Node: 158 + -2.558688446134715e-02 + -1.841019731224498e-01 + 4.638007494765678e-01 +Node: 159 + -2.468222564382960e-03 + -5.100631182990298e-02 + 4.970701094227351e-01 +Node: 160 + 4.467609383512137e-02 + -1.608178737905387e-01 + 4.713730776792378e-01 +Node: 161 + 6.383793283641934e-02 + -1.431724856883514e-01 + 4.743551281852820e-01 +Node: 162 + 1.163393476697186e-01 + -1.711938748360189e-01 + 4.549287173574261e-01 +#Fields=2 +1) marker_data_coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=1 (value) + y. #Values=1 (value) + z. #Values=1 (value) +2) marker_data_name, field, rectangular cartesian, string, #Components=1 + 1. #Values=1 (value) +Node: 1000001 + 0.000000000000000e+00 + 0.000000000000000e+00 + -5.000000000000000e-01 + bottom +Node: 1000002 + 0.000000000000000e+00 + 0.000000000000000e+00 + 5.000000000000000e-01 + top +Node: 1000003 + -5.000000000000000e-01 + 0.000000000000000e+00 + 0.000000000000000e+00 + left +Node: 1000004 + 5.000000000000000e-01 + 0.000000000000000e+00 + 0.000000000000000e+00 + right + Group name: bottom +!#nodeset datapoints +Shape. Dimension=0 +#Fields=0 +Node: 1 +Node: 2 +Node: 3 +Node: 4 +Node: 5 +Node: 6 +Node: 7 +Node: 8 +Node: 9 +Node: 10 +Node: 11 +Node: 13 +Node: 14 +Node: 17 +Node: 18 +Node: 19 +Node: 21 +Node: 22 +Node: 24 +Node: 25 +Node: 26 +Node: 28 +Node: 30 +Node: 32 +Node: 34 + Group name: marker +!#nodeset datapoints +Shape. Dimension=0 +#Fields=0 +Node: 1000001 +Node: 1000002 +Node: 1000003 +Node: 1000004 + Group name: sides +!#nodeset datapoints +Shape. Dimension=0 +#Fields=0 +Node: 12 +Node: 15 +Node: 16 +Node: 20 +Node: 23 +Node: 27 +Node: 29 +Node: 31 +Node: 33 +Node: 35 +Node: 36 +Node: 37 +Node: 38 +Node: 39 +Node: 40 +Node: 41 +Node: 42 +Node: 43 +Node: 44 +Node: 45 +Node: 46 +Node: 47 +Node: 48 +Node: 49 +Node: 50 +Node: 51 +Node: 52 +Node: 53 +Node: 54 +Node: 55 +Node: 56 +Node: 57 +Node: 58 +Node: 59 +Node: 60 +Node: 61 +Node: 62 +Node: 63 +Node: 64 +Node: 65 +Node: 66 +Node: 67 +Node: 68 +Node: 69 +Node: 70 +Node: 71 +Node: 72 +Node: 73 +Node: 74 +Node: 75 +Node: 76 +Node: 77 +Node: 78 +Node: 79 +Node: 80 +Node: 81 +Node: 82 +Node: 83 +Node: 84 +Node: 85 +Node: 86 +Node: 87 +Node: 88 +Node: 89 +Node: 90 +Node: 91 +Node: 92 +Node: 93 +Node: 94 +Node: 95 +Node: 96 +Node: 97 +Node: 98 +Node: 99 +Node: 100 +Node: 101 +Node: 102 +Node: 103 +Node: 104 +Node: 105 +Node: 106 +Node: 107 +Node: 108 +Node: 109 +Node: 110 +Node: 111 +Node: 112 +Node: 113 +Node: 114 +Node: 115 +Node: 116 +Node: 117 +Node: 118 +Node: 119 +Node: 123 +Node: 124 +Node: 126 +Node: 130 +Node: 131 +Node: 133 +Node: 134 +Node: 139 +Node: 140 +Node: 141 +Node: 143 +Node: 144 +Node: 145 +Node: 146 +Node: 147 +Node: 148 +Node: 149 + Group name: top +!#nodeset datapoints +Shape. Dimension=0 +#Fields=0 +Node: 120 +Node: 121 +Node: 122 +Node: 125 +Node: 127 +Node: 128 +Node: 129 +Node: 132 +Node: 135 +Node: 136 +Node: 137 +Node: 138 +Node: 142 +Node: 150 +Node: 151 +Node: 152 +Node: 153 +Node: 154 +Node: 155 +Node: 156 +Node: 157 +Node: 158 +Node: 159 +Node: 160 +Node: 161 +Node: 162 diff --git a/tests/resources/cube_to_sphere_data_regular.exf b/tests/resources/cube_to_sphere_data_regular.exf new file mode 100644 index 0000000..ec388a7 --- /dev/null +++ b/tests/resources/cube_to_sphere_data_regular.exf @@ -0,0 +1,1497 @@ +EX Version: 2 +Region: / +!#nodeset nodes +Shape. Dimension=0 +#Fields=1 +1) data_coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=1 (value) + y. #Values=1 (value) + z. #Values=1 (value) +Node: 1 + 6.481555104255676e-02 + 9.071839042007923e-03 + -4.956691563129425e-01 +Node: 2 + 6.066783517599106e-02 + 2.455133199691772e-02 + -4.956691563129425e-01 +Node: 3 + 1.900602132081985e-01 + 2.663808315992355e-02 + -4.618872404098511e-01 +Node: 4 + 1.779160201549530e-01 + 7.196085155010223e-02 + -4.618872404098511e-01 +Node: 5 + 5.159599706530571e-02 + 4.026421904563904e-02 + -4.956691563129425e-01 +Node: 6 + 4.026421904563904e-02 + 5.159599706530571e-02 + -4.956691563129425e-01 +Node: 7 + 1.512779295444489e-01 + 1.180993616580963e-01 + -4.618872404098511e-01 +Node: 8 + 1.180993616580963e-01 + 1.512779295444489e-01 + -4.618872404098511e-01 +Node: 9 + 2.455133199691772e-02 + 6.066783517599106e-02 + -4.956691563129425e-01 +Node: 10 + 9.071839042007923e-03 + 6.481555104255676e-02 + -4.956691563129425e-01 +Node: 11 + 7.196085155010223e-02 + 1.779160201549530e-01 + -4.618872404098511e-01 +Node: 12 + 2.663808315992355e-02 + 1.900602132081985e-01 + -4.618872404098511e-01 +Node: 13 + -9.071839042007923e-03 + 6.481555104255676e-02 + -4.956691563129425e-01 +Node: 14 + -2.455133199691772e-02 + 6.066783517599106e-02 + -4.956691563129425e-01 +Node: 15 + -2.663808315992355e-02 + 1.900602132081985e-01 + -4.618872404098511e-01 +Node: 16 + -7.196085155010223e-02 + 1.779160201549530e-01 + -4.618872404098511e-01 +Node: 17 + -4.026421904563904e-02 + 5.159599706530571e-02 + -4.956691563129425e-01 +Node: 18 + -5.159599706530571e-02 + 4.026421904563904e-02 + -4.956691563129425e-01 +Node: 19 + -1.180993616580963e-01 + 1.512779295444489e-01 + -4.618872404098511e-01 +Node: 20 + -1.512779295444489e-01 + 1.180993616580963e-01 + -4.618872404098511e-01 +Node: 21 + -6.066783517599106e-02 + 2.455133199691772e-02 + -4.956691563129425e-01 +Node: 22 + -6.481555104255676e-02 + 9.071839042007923e-03 + -4.956691563129425e-01 +Node: 23 + -1.779160201549530e-01 + 7.196085155010223e-02 + -4.618872404098511e-01 +Node: 24 + -1.900602132081985e-01 + 2.663808315992355e-02 + -4.618872404098511e-01 +Node: 25 + -6.481555104255676e-02 + -9.071839042007923e-03 + -4.956691563129425e-01 +Node: 26 + -6.066783517599106e-02 + -2.455133199691772e-02 + -4.956691563129425e-01 +Node: 27 + -1.900602132081985e-01 + -2.663808315992355e-02 + -4.618872404098511e-01 +Node: 28 + -1.779160201549530e-01 + -7.196085155010223e-02 + -4.618872404098511e-01 +Node: 29 + -5.159599706530571e-02 + -4.026421904563904e-02 + -4.956691563129425e-01 +Node: 30 + -4.026421904563904e-02 + -5.159599706530571e-02 + -4.956691563129425e-01 +Node: 31 + -1.512779295444489e-01 + -1.180993616580963e-01 + -4.618872404098511e-01 +Node: 32 + -1.180993616580963e-01 + -1.512779295444489e-01 + -4.618872404098511e-01 +Node: 33 + -2.455133199691772e-02 + -6.066783517599106e-02 + -4.956691563129425e-01 +Node: 34 + -9.071839042007923e-03 + -6.481555104255676e-02 + -4.956691563129425e-01 +Node: 35 + -7.196085155010223e-02 + -1.779160201549530e-01 + -4.618872404098511e-01 +Node: 36 + -2.663808315992355e-02 + -1.900602132081985e-01 + -4.618872404098511e-01 +Node: 37 + 9.071839042007923e-03 + -6.481555104255676e-02 + -4.956691563129425e-01 +Node: 38 + 2.455133199691772e-02 + -6.066783517599106e-02 + -4.956691563129425e-01 +Node: 39 + 2.663808315992355e-02 + -1.900602132081985e-01 + -4.618872404098511e-01 +Node: 40 + 7.196085155010223e-02 + -1.779160201549530e-01 + -4.618872404098511e-01 +Node: 41 + 4.026421904563904e-02 + -5.159599706530571e-02 + -4.956691563129425e-01 +Node: 42 + 5.159599706530571e-02 + -4.026421904563904e-02 + -4.956691563129425e-01 +Node: 43 + 1.180993616580963e-01 + -1.512779295444489e-01 + -4.618872404098511e-01 +Node: 44 + 1.512779295444489e-01 + -1.180993616580963e-01 + -4.618872404098511e-01 +Node: 45 + 6.066783517599106e-02 + -2.455133199691772e-02 + -4.956691563129425e-01 +Node: 46 + 6.481555104255676e-02 + -9.071839042007923e-03 + -4.956691563129425e-01 +Node: 47 + 1.779160201549530e-01 + -7.196085155010223e-02 + -4.618872404098511e-01 +Node: 48 + 1.900602132081985e-01 + -2.663808315992355e-02 + -4.618872404098511e-01 +Node: 49 + 3.013909161090851e-01 + 3.836841881275177e-02 + -3.966369330883026e-01 +Node: 50 + 2.801963984966278e-01 + 1.174674332141876e-01 + -3.966369330883026e-01 +Node: 51 + 3.932968676090240e-01 + 5.216884985566139e-02 + -3.043430149555206e-01 +Node: 52 + 3.666895031929016e-01 + 1.514688879251480e-01 + -3.043430149555206e-01 +Node: 53 + 2.418279796838760e-01 + 1.839234828948975e-01 + -3.966369330883026e-01 +Node: 54 + 1.839234828948975e-01 + 2.418279796838760e-01 + -3.966369330883026e-01 +Node: 55 + 3.145206570625305e-01 + 2.418279796838760e-01 + -3.043430149555206e-01 +Node: 56 + 2.418279796838760e-01 + 3.145206570625305e-01 + -3.043430149555206e-01 +Node: 57 + 1.174674332141876e-01 + 2.801963984966278e-01 + -3.966369330883026e-01 +Node: 58 + 3.836841881275177e-02 + 3.013909161090851e-01 + -3.966369330883026e-01 +Node: 59 + 1.514688879251480e-01 + 3.666895031929016e-01 + -3.043430149555206e-01 +Node: 60 + 5.216884985566139e-02 + 3.932968676090240e-01 + -3.043430149555206e-01 +Node: 61 + -3.836841881275177e-02 + 3.013909161090851e-01 + -3.966369330883026e-01 +Node: 62 + -1.174674332141876e-01 + 2.801963984966278e-01 + -3.966369330883026e-01 +Node: 63 + -5.216884985566139e-02 + 3.932968676090240e-01 + -3.043430149555206e-01 +Node: 64 + -1.514688879251480e-01 + 3.666895031929016e-01 + -3.043430149555206e-01 +Node: 65 + -1.839234828948975e-01 + 2.418279796838760e-01 + -3.966369330883026e-01 +Node: 66 + -2.418279796838760e-01 + 1.839234828948975e-01 + -3.966369330883026e-01 +Node: 67 + -2.418279796838760e-01 + 3.145206570625305e-01 + -3.043430149555206e-01 +Node: 68 + -3.145206570625305e-01 + 2.418279796838760e-01 + -3.043430149555206e-01 +Node: 69 + -2.801963984966278e-01 + 1.174674332141876e-01 + -3.966369330883026e-01 +Node: 70 + -3.013909161090851e-01 + 3.836841881275177e-02 + -3.966369330883026e-01 +Node: 71 + -3.666895031929016e-01 + 1.514688879251480e-01 + -3.043430149555206e-01 +Node: 72 + -3.932968676090240e-01 + 5.216884985566139e-02 + -3.043430149555206e-01 +Node: 73 + -3.013909161090851e-01 + -3.836841881275177e-02 + -3.966369330883026e-01 +Node: 74 + -2.801963984966278e-01 + -1.174674332141876e-01 + -3.966369330883026e-01 +Node: 75 + -3.932968676090240e-01 + -5.216884985566139e-02 + -3.043430149555206e-01 +Node: 76 + -3.666895031929016e-01 + -1.514688879251480e-01 + -3.043430149555206e-01 +Node: 77 + -2.418279796838760e-01 + -1.839234828948975e-01 + -3.966369330883026e-01 +Node: 78 + -1.839234828948975e-01 + -2.418279796838760e-01 + -3.966369330883026e-01 +Node: 79 + -3.145206570625305e-01 + -2.418279796838760e-01 + -3.043430149555206e-01 +Node: 80 + -2.418279796838760e-01 + -3.145206570625305e-01 + -3.043430149555206e-01 +Node: 81 + -1.174674332141876e-01 + -2.801963984966278e-01 + -3.966369330883026e-01 +Node: 82 + -3.836841881275177e-02 + -3.013909161090851e-01 + -3.966369330883026e-01 +Node: 83 + -1.514688879251480e-01 + -3.666895031929016e-01 + -3.043430149555206e-01 +Node: 84 + -5.216884985566139e-02 + -3.932968676090240e-01 + -3.043430149555206e-01 +Node: 85 + 3.836841881275177e-02 + -3.013909161090851e-01 + -3.966369330883026e-01 +Node: 86 + 1.174674332141876e-01 + -2.801963984966278e-01 + -3.966369330883026e-01 +Node: 87 + 5.216884985566139e-02 + -3.932968676090240e-01 + -3.043430149555206e-01 +Node: 88 + 1.514688879251480e-01 + -3.666895031929016e-01 + -3.043430149555206e-01 +Node: 89 + 1.839234828948975e-01 + -2.418279796838760e-01 + -3.966369330883026e-01 +Node: 90 + 2.418279796838760e-01 + -1.839234828948975e-01 + -3.966369330883026e-01 +Node: 91 + 2.418279796838760e-01 + -3.145206570625305e-01 + -3.043430149555206e-01 +Node: 92 + 3.145206570625305e-01 + -2.418279796838760e-01 + -3.043430149555206e-01 +Node: 93 + 2.801963984966278e-01 + -1.174674332141876e-01 + -3.966369330883026e-01 +Node: 94 + 3.013909161090851e-01 + -3.836841881275177e-02 + -3.966369330883026e-01 +Node: 95 + 3.666895031929016e-01 + -1.514688879251480e-01 + -3.043430149555206e-01 +Node: 96 + 3.932968676090240e-01 + -5.216884985566139e-02 + -3.043430149555206e-01 +Node: 97 + 4.576606154441833e-01 + 5.931245163083076e-02 + -1.913261562585831e-01 +Node: 98 + 4.260019361972809e-01 + 1.774642169475555e-01 + -1.913261562585831e-01 +Node: 99 + 4.913005232810974e-01 + 6.436375528573990e-02 + -6.525030732154846e-02 +Node: 100 + 4.576606154441833e-01 + 1.899096220731735e-01 + -6.525030732154846e-02 +Node: 101 + 3.666895031929016e-01 + 2.801963984966278e-01 + -1.913261562585831e-01 +Node: 102 + 2.801963984966278e-01 + 3.666895031929016e-01 + -1.913261562585831e-01 +Node: 103 + 3.932968676090240e-01 + 3.013909161090851e-01 + -6.525030732154846e-02 +Node: 104 + 3.013909161090851e-01 + 3.932968676090240e-01 + -6.525030732154846e-02 +Node: 105 + 1.774642169475555e-01 + 4.260019361972809e-01 + -1.913261562585831e-01 +Node: 106 + 5.931245163083076e-02 + 4.576606154441833e-01 + -1.913261562585831e-01 +Node: 107 + 1.899096220731735e-01 + 4.576606154441833e-01 + -6.525030732154846e-02 +Node: 108 + 6.436375528573990e-02 + 4.913005232810974e-01 + -6.525030732154846e-02 +Node: 109 + -5.931245163083076e-02 + 4.576606154441833e-01 + -1.913261562585831e-01 +Node: 110 + -1.774642169475555e-01 + 4.260019361972809e-01 + -1.913261562585831e-01 +Node: 111 + -6.436375528573990e-02 + 4.913005232810974e-01 + -6.525030732154846e-02 +Node: 112 + -1.899096220731735e-01 + 4.576606154441833e-01 + -6.525030732154846e-02 +Node: 113 + -2.801963984966278e-01 + 3.666895031929016e-01 + -1.913261562585831e-01 +Node: 114 + -3.666895031929016e-01 + 2.801963984966278e-01 + -1.913261562585831e-01 +Node: 115 + -3.013909161090851e-01 + 3.932968676090240e-01 + -6.525030732154846e-02 +Node: 116 + -3.932968676090240e-01 + 3.013909161090851e-01 + -6.525030732154846e-02 +Node: 117 + -4.260019361972809e-01 + 1.774642169475555e-01 + -1.913261562585831e-01 +Node: 118 + -4.576606154441833e-01 + 5.931245163083076e-02 + -1.913261562585831e-01 +Node: 119 + -4.576606154441833e-01 + 1.899096220731735e-01 + -6.525030732154846e-02 +Node: 120 + -4.913005232810974e-01 + 6.436375528573990e-02 + -6.525030732154846e-02 +Node: 121 + -4.576606154441833e-01 + -5.931245163083076e-02 + -1.913261562585831e-01 +Node: 122 + -4.260019361972809e-01 + -1.774642169475555e-01 + -1.913261562585831e-01 +Node: 123 + -4.913005232810974e-01 + -6.436375528573990e-02 + -6.525030732154846e-02 +Node: 124 + -4.576606154441833e-01 + -1.899096220731735e-01 + -6.525030732154846e-02 +Node: 125 + -3.666895031929016e-01 + -2.801963984966278e-01 + -1.913261562585831e-01 +Node: 126 + -2.801963984966278e-01 + -3.666895031929016e-01 + -1.913261562585831e-01 +Node: 127 + -3.932968676090240e-01 + -3.013909161090851e-01 + -6.525030732154846e-02 +Node: 128 + -3.013909161090851e-01 + -3.932968676090240e-01 + -6.525030732154846e-02 +Node: 129 + -1.774642169475555e-01 + -4.260019361972809e-01 + -1.913261562585831e-01 +Node: 130 + -5.931245163083076e-02 + -4.576606154441833e-01 + -1.913261562585831e-01 +Node: 131 + -1.899096220731735e-01 + -4.576606154441833e-01 + -6.525030732154846e-02 +Node: 132 + -6.436375528573990e-02 + -4.913005232810974e-01 + -6.525030732154846e-02 +Node: 133 + 5.931245163083076e-02 + -4.576606154441833e-01 + -1.913261562585831e-01 +Node: 134 + 1.774642169475555e-01 + -4.260019361972809e-01 + -1.913261562585831e-01 +Node: 135 + 6.436375528573990e-02 + -4.913005232810974e-01 + -6.525030732154846e-02 +Node: 136 + 1.899096220731735e-01 + -4.576606154441833e-01 + -6.525030732154846e-02 +Node: 137 + 2.801963984966278e-01 + -3.666895031929016e-01 + -1.913261562585831e-01 +Node: 138 + 3.666895031929016e-01 + -2.801963984966278e-01 + -1.913261562585831e-01 +Node: 139 + 3.013909161090851e-01 + -3.932968676090240e-01 + -6.525030732154846e-02 +Node: 140 + 3.932968676090240e-01 + -3.013909161090851e-01 + -6.525030732154846e-02 +Node: 141 + 4.260019361972809e-01 + -1.774642169475555e-01 + -1.913261562585831e-01 +Node: 142 + 4.576606154441833e-01 + -5.931245163083076e-02 + -1.913261562585831e-01 +Node: 143 + 4.576606154441833e-01 + -1.899096220731735e-01 + -6.525030732154846e-02 +Node: 144 + 4.913005232810974e-01 + -6.436375528573990e-02 + -6.525030732154846e-02 +Node: 145 + 4.913005232810974e-01 + 6.436375528573990e-02 + 6.525030732154846e-02 +Node: 146 + 4.576606154441833e-01 + 1.899096220731735e-01 + 6.525030732154846e-02 +Node: 147 + 4.576606154441833e-01 + 5.931245163083076e-02 + 1.913261562585831e-01 +Node: 148 + 4.260019361972809e-01 + 1.774642169475555e-01 + 1.913261562585831e-01 +Node: 149 + 3.932968676090240e-01 + 3.013909161090851e-01 + 6.525030732154846e-02 +Node: 150 + 3.013909161090851e-01 + 3.932968676090240e-01 + 6.525030732154846e-02 +Node: 151 + 3.666895031929016e-01 + 2.801963984966278e-01 + 1.913261562585831e-01 +Node: 152 + 2.801963984966278e-01 + 3.666895031929016e-01 + 1.913261562585831e-01 +Node: 153 + 1.899096220731735e-01 + 4.576606154441833e-01 + 6.525030732154846e-02 +Node: 154 + 6.436375528573990e-02 + 4.913005232810974e-01 + 6.525030732154846e-02 +Node: 155 + 1.774642169475555e-01 + 4.260019361972809e-01 + 1.913261562585831e-01 +Node: 156 + 5.931245163083076e-02 + 4.576606154441833e-01 + 1.913261562585831e-01 +Node: 157 + -6.436375528573990e-02 + 4.913005232810974e-01 + 6.525030732154846e-02 +Node: 158 + -1.899096220731735e-01 + 4.576606154441833e-01 + 6.525030732154846e-02 +Node: 159 + -5.931245163083076e-02 + 4.576606154441833e-01 + 1.913261562585831e-01 +Node: 160 + -1.774642169475555e-01 + 4.260019361972809e-01 + 1.913261562585831e-01 +Node: 161 + -3.013909161090851e-01 + 3.932968676090240e-01 + 6.525030732154846e-02 +Node: 162 + -3.932968676090240e-01 + 3.013909161090851e-01 + 6.525030732154846e-02 +Node: 163 + -2.801963984966278e-01 + 3.666895031929016e-01 + 1.913261562585831e-01 +Node: 164 + -3.666895031929016e-01 + 2.801963984966278e-01 + 1.913261562585831e-01 +Node: 165 + -4.576606154441833e-01 + 1.899096220731735e-01 + 6.525030732154846e-02 +Node: 166 + -4.913005232810974e-01 + 6.436375528573990e-02 + 6.525030732154846e-02 +Node: 167 + -4.260019361972809e-01 + 1.774642169475555e-01 + 1.913261562585831e-01 +Node: 168 + -4.576606154441833e-01 + 5.931245163083076e-02 + 1.913261562585831e-01 +Node: 169 + -4.913005232810974e-01 + -6.436375528573990e-02 + 6.525030732154846e-02 +Node: 170 + -4.576606154441833e-01 + -1.899096220731735e-01 + 6.525030732154846e-02 +Node: 171 + -4.576606154441833e-01 + -5.931245163083076e-02 + 1.913261562585831e-01 +Node: 172 + -4.260019361972809e-01 + -1.774642169475555e-01 + 1.913261562585831e-01 +Node: 173 + -3.932968676090240e-01 + -3.013909161090851e-01 + 6.525030732154846e-02 +Node: 174 + -3.013909161090851e-01 + -3.932968676090240e-01 + 6.525030732154846e-02 +Node: 175 + -3.666895031929016e-01 + -2.801963984966278e-01 + 1.913261562585831e-01 +Node: 176 + -2.801963984966278e-01 + -3.666895031929016e-01 + 1.913261562585831e-01 +Node: 177 + -1.899096220731735e-01 + -4.576606154441833e-01 + 6.525030732154846e-02 +Node: 178 + -6.436375528573990e-02 + -4.913005232810974e-01 + 6.525030732154846e-02 +Node: 179 + -1.774642169475555e-01 + -4.260019361972809e-01 + 1.913261562585831e-01 +Node: 180 + -5.931245163083076e-02 + -4.576606154441833e-01 + 1.913261562585831e-01 +Node: 181 + 6.436375528573990e-02 + -4.913005232810974e-01 + 6.525030732154846e-02 +Node: 182 + 1.899096220731735e-01 + -4.576606154441833e-01 + 6.525030732154846e-02 +Node: 183 + 5.931245163083076e-02 + -4.576606154441833e-01 + 1.913261562585831e-01 +Node: 184 + 1.774642169475555e-01 + -4.260019361972809e-01 + 1.913261562585831e-01 +Node: 185 + 3.013909161090851e-01 + -3.932968676090240e-01 + 6.525030732154846e-02 +Node: 186 + 3.932968676090240e-01 + -3.013909161090851e-01 + 6.525030732154846e-02 +Node: 187 + 2.801963984966278e-01 + -3.666895031929016e-01 + 1.913261562585831e-01 +Node: 188 + 3.666895031929016e-01 + -2.801963984966278e-01 + 1.913261562585831e-01 +Node: 189 + 4.576606154441833e-01 + -1.899096220731735e-01 + 6.525030732154846e-02 +Node: 190 + 4.913005232810974e-01 + -6.436375528573990e-02 + 6.525030732154846e-02 +Node: 191 + 4.260019361972809e-01 + -1.774642169475555e-01 + 1.913261562585831e-01 +Node: 192 + 4.576606154441833e-01 + -5.931245163083076e-02 + 1.913261562585831e-01 +Node: 193 + 3.932968676090240e-01 + 5.216884985566139e-02 + 3.043430149555206e-01 +Node: 194 + 3.666895031929016e-01 + 1.514688879251480e-01 + 3.043430149555206e-01 +Node: 195 + 3.013909161090851e-01 + 3.836841881275177e-02 + 3.966369330883026e-01 +Node: 196 + 2.801963984966278e-01 + 1.174674332141876e-01 + 3.966369330883026e-01 +Node: 197 + 3.145206570625305e-01 + 2.418279796838760e-01 + 3.043430149555206e-01 +Node: 198 + 2.418279796838760e-01 + 3.145206570625305e-01 + 3.043430149555206e-01 +Node: 199 + 2.418279796838760e-01 + 1.839234828948975e-01 + 3.966369330883026e-01 +Node: 200 + 1.839234828948975e-01 + 2.418279796838760e-01 + 3.966369330883026e-01 +Node: 201 + 1.514688879251480e-01 + 3.666895031929016e-01 + 3.043430149555206e-01 +Node: 202 + 5.216884985566139e-02 + 3.932968676090240e-01 + 3.043430149555206e-01 +Node: 203 + 1.174674332141876e-01 + 2.801963984966278e-01 + 3.966369330883026e-01 +Node: 204 + 3.836841881275177e-02 + 3.013909161090851e-01 + 3.966369330883026e-01 +Node: 205 + -5.216884985566139e-02 + 3.932968676090240e-01 + 3.043430149555206e-01 +Node: 206 + -1.514688879251480e-01 + 3.666895031929016e-01 + 3.043430149555206e-01 +Node: 207 + -3.836841881275177e-02 + 3.013909161090851e-01 + 3.966369330883026e-01 +Node: 208 + -1.174674332141876e-01 + 2.801963984966278e-01 + 3.966369330883026e-01 +Node: 209 + -2.418279796838760e-01 + 3.145206570625305e-01 + 3.043430149555206e-01 +Node: 210 + -3.145206570625305e-01 + 2.418279796838760e-01 + 3.043430149555206e-01 +Node: 211 + -1.839234828948975e-01 + 2.418279796838760e-01 + 3.966369330883026e-01 +Node: 212 + -2.418279796838760e-01 + 1.839234828948975e-01 + 3.966369330883026e-01 +Node: 213 + -3.666895031929016e-01 + 1.514688879251480e-01 + 3.043430149555206e-01 +Node: 214 + -3.932968676090240e-01 + 5.216884985566139e-02 + 3.043430149555206e-01 +Node: 215 + -2.801963984966278e-01 + 1.174674332141876e-01 + 3.966369330883026e-01 +Node: 216 + -3.013909161090851e-01 + 3.836841881275177e-02 + 3.966369330883026e-01 +Node: 217 + -3.932968676090240e-01 + -5.216884985566139e-02 + 3.043430149555206e-01 +Node: 218 + -3.666895031929016e-01 + -1.514688879251480e-01 + 3.043430149555206e-01 +Node: 219 + -3.013909161090851e-01 + -3.836841881275177e-02 + 3.966369330883026e-01 +Node: 220 + -2.801963984966278e-01 + -1.174674332141876e-01 + 3.966369330883026e-01 +Node: 221 + -3.145206570625305e-01 + -2.418279796838760e-01 + 3.043430149555206e-01 +Node: 222 + -2.418279796838760e-01 + -3.145206570625305e-01 + 3.043430149555206e-01 +Node: 223 + -2.418279796838760e-01 + -1.839234828948975e-01 + 3.966369330883026e-01 +Node: 224 + -1.839234828948975e-01 + -2.418279796838760e-01 + 3.966369330883026e-01 +Node: 225 + -1.514688879251480e-01 + -3.666895031929016e-01 + 3.043430149555206e-01 +Node: 226 + -5.216884985566139e-02 + -3.932968676090240e-01 + 3.043430149555206e-01 +Node: 227 + -1.174674332141876e-01 + -2.801963984966278e-01 + 3.966369330883026e-01 +Node: 228 + -3.836841881275177e-02 + -3.013909161090851e-01 + 3.966369330883026e-01 +Node: 229 + 5.216884985566139e-02 + -3.932968676090240e-01 + 3.043430149555206e-01 +Node: 230 + 1.514688879251480e-01 + -3.666895031929016e-01 + 3.043430149555206e-01 +Node: 231 + 3.836841881275177e-02 + -3.013909161090851e-01 + 3.966369330883026e-01 +Node: 232 + 1.174674332141876e-01 + -2.801963984966278e-01 + 3.966369330883026e-01 +Node: 233 + 2.418279796838760e-01 + -3.145206570625305e-01 + 3.043430149555206e-01 +Node: 234 + 3.145206570625305e-01 + -2.418279796838760e-01 + 3.043430149555206e-01 +Node: 235 + 1.839234828948975e-01 + -2.418279796838760e-01 + 3.966369330883026e-01 +Node: 236 + 2.418279796838760e-01 + -1.839234828948975e-01 + 3.966369330883026e-01 +Node: 237 + 3.666895031929016e-01 + -1.514688879251480e-01 + 3.043430149555206e-01 +Node: 238 + 3.932968676090240e-01 + -5.216884985566139e-02 + 3.043430149555206e-01 +Node: 239 + 2.801963984966278e-01 + -1.174674332141876e-01 + 3.966369330883026e-01 +Node: 240 + 3.013909161090851e-01 + -3.836841881275177e-02 + 3.966369330883026e-01 +Node: 241 + 1.900602132081985e-01 + 2.663808315992355e-02 + 4.618872404098511e-01 +Node: 242 + 1.779160201549530e-01 + 7.196085155010223e-02 + 4.618872404098511e-01 +Node: 243 + 6.481555104255676e-02 + 9.071839042007923e-03 + 4.956691563129425e-01 +Node: 244 + 6.066783517599106e-02 + 2.455133199691772e-02 + 4.956691563129425e-01 +Node: 245 + 1.512779295444489e-01 + 1.180993616580963e-01 + 4.618872404098511e-01 +Node: 246 + 1.180993616580963e-01 + 1.512779295444489e-01 + 4.618872404098511e-01 +Node: 247 + 5.159599706530571e-02 + 4.026421904563904e-02 + 4.956691563129425e-01 +Node: 248 + 4.026421904563904e-02 + 5.159599706530571e-02 + 4.956691563129425e-01 +Node: 249 + 7.196085155010223e-02 + 1.779160201549530e-01 + 4.618872404098511e-01 +Node: 250 + 2.663808315992355e-02 + 1.900602132081985e-01 + 4.618872404098511e-01 +Node: 251 + 2.455133199691772e-02 + 6.066783517599106e-02 + 4.956691563129425e-01 +Node: 252 + 9.071839042007923e-03 + 6.481555104255676e-02 + 4.956691563129425e-01 +Node: 253 + -2.663808315992355e-02 + 1.900602132081985e-01 + 4.618872404098511e-01 +Node: 254 + -7.196085155010223e-02 + 1.779160201549530e-01 + 4.618872404098511e-01 +Node: 255 + -9.071839042007923e-03 + 6.481555104255676e-02 + 4.956691563129425e-01 +Node: 256 + -2.455133199691772e-02 + 6.066783517599106e-02 + 4.956691563129425e-01 +Node: 257 + -1.180993616580963e-01 + 1.512779295444489e-01 + 4.618872404098511e-01 +Node: 258 + -1.512779295444489e-01 + 1.180993616580963e-01 + 4.618872404098511e-01 +Node: 259 + -4.026421904563904e-02 + 5.159599706530571e-02 + 4.956691563129425e-01 +Node: 260 + -5.159599706530571e-02 + 4.026421904563904e-02 + 4.956691563129425e-01 +Node: 261 + -1.779160201549530e-01 + 7.196085155010223e-02 + 4.618872404098511e-01 +Node: 262 + -1.900602132081985e-01 + 2.663808315992355e-02 + 4.618872404098511e-01 +Node: 263 + -6.066783517599106e-02 + 2.455133199691772e-02 + 4.956691563129425e-01 +Node: 264 + -6.481555104255676e-02 + 9.071839042007923e-03 + 4.956691563129425e-01 +Node: 265 + -1.900602132081985e-01 + -2.663808315992355e-02 + 4.618872404098511e-01 +Node: 266 + -1.779160201549530e-01 + -7.196085155010223e-02 + 4.618872404098511e-01 +Node: 267 + -6.481555104255676e-02 + -9.071839042007923e-03 + 4.956691563129425e-01 +Node: 268 + -6.066783517599106e-02 + -2.455133199691772e-02 + 4.956691563129425e-01 +Node: 269 + -1.512779295444489e-01 + -1.180993616580963e-01 + 4.618872404098511e-01 +Node: 270 + -1.180993616580963e-01 + -1.512779295444489e-01 + 4.618872404098511e-01 +Node: 271 + -5.159599706530571e-02 + -4.026421904563904e-02 + 4.956691563129425e-01 +Node: 272 + -4.026421904563904e-02 + -5.159599706530571e-02 + 4.956691563129425e-01 +Node: 273 + -7.196085155010223e-02 + -1.779160201549530e-01 + 4.618872404098511e-01 +Node: 274 + -2.663808315992355e-02 + -1.900602132081985e-01 + 4.618872404098511e-01 +Node: 275 + -2.455133199691772e-02 + -6.066783517599106e-02 + 4.956691563129425e-01 +Node: 276 + -9.071839042007923e-03 + -6.481555104255676e-02 + 4.956691563129425e-01 +Node: 277 + 2.663808315992355e-02 + -1.900602132081985e-01 + 4.618872404098511e-01 +Node: 278 + 7.196085155010223e-02 + -1.779160201549530e-01 + 4.618872404098511e-01 +Node: 279 + 9.071839042007923e-03 + -6.481555104255676e-02 + 4.956691563129425e-01 +Node: 280 + 2.455133199691772e-02 + -6.066783517599106e-02 + 4.956691563129425e-01 +Node: 281 + 1.180993616580963e-01 + -1.512779295444489e-01 + 4.618872404098511e-01 +Node: 282 + 1.512779295444489e-01 + -1.180993616580963e-01 + 4.618872404098511e-01 +Node: 283 + 4.026421904563904e-02 + -5.159599706530571e-02 + 4.956691563129425e-01 +Node: 284 + 5.159599706530571e-02 + -4.026421904563904e-02 + 4.956691563129425e-01 +Node: 285 + 1.779160201549530e-01 + -7.196085155010223e-02 + 4.618872404098511e-01 +Node: 286 + 1.900602132081985e-01 + -2.663808315992355e-02 + 4.618872404098511e-01 +Node: 287 + 6.066783517599106e-02 + -2.455133199691772e-02 + 4.956691563129425e-01 +Node: 288 + 6.481555104255676e-02 + -9.071839042007923e-03 + 4.956691563129425e-01 +!#nodeset datapoints +#Fields=2 +1) marker_data_coordinates, coordinate, rectangular cartesian, real, #Components=3 + x. #Values=1 (value) + y. #Values=1 (value) + z. #Values=1 (value) +2) marker_data_name, field, rectangular cartesian, string, #Components=1 + 1. #Values=1 (value) +Node: 1 + 0.000000000000000e+00 + 0.000000000000000e+00 + -5.000000000000000e-01 + bottom +Node: 2 + 0.000000000000000e+00 + 0.000000000000000e+00 + 5.000000000000000e-01 + top +Node: 3 + -5.000000000000000e-01 + 0.000000000000000e+00 + 0.000000000000000e+00 + left +Node: 4 + 5.000000000000000e-01 + 0.000000000000000e+00 + 0.000000000000000e+00 + right + Group name: bottom +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 1 +Node: 2 +Node: 3 +Node: 4 +Node: 5 +Node: 6 +Node: 7 +Node: 8 +Node: 9 +Node: 10 +Node: 11 +Node: 12 +Node: 13 +Node: 14 +Node: 15 +Node: 16 +Node: 17 +Node: 18 +Node: 19 +Node: 20 +Node: 21 +Node: 22 +Node: 23 +Node: 24 +Node: 25 +Node: 26 +Node: 27 +Node: 28 +Node: 29 +Node: 30 +Node: 31 +Node: 32 +Node: 33 +Node: 34 +Node: 35 +Node: 36 +Node: 37 +Node: 38 +Node: 39 +Node: 40 +Node: 41 +Node: 42 +Node: 43 +Node: 44 +Node: 45 +Node: 46 +Node: 47 +Node: 48 +Node: 49 +Node: 50 +Node: 53 +Node: 54 +Node: 57 +Node: 58 +Node: 61 +Node: 62 +Node: 65 +Node: 66 +Node: 69 +Node: 70 +Node: 73 +Node: 74 +Node: 77 +Node: 78 +Node: 81 +Node: 82 +Node: 85 +Node: 86 +Node: 89 +Node: 90 +Node: 93 +Node: 94 + Group name: marker +!#nodeset datapoints +Shape. Dimension=0 +#Fields=0 +Node: 1 +Node: 2 +Node: 3 +Node: 4 + Group name: sides +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 51 +Node: 52 +Node: 55 +Node: 56 +Node: 59 +Node: 60 +Node: 63 +Node: 64 +Node: 67 +Node: 68 +Node: 71 +Node: 72 +Node: 75 +Node: 76 +Node: 79 +Node: 80 +Node: 83 +Node: 84 +Node: 87 +Node: 88 +Node: 91 +Node: 92 +Node: 95 +Node: 96 +Node: 97 +Node: 98 +Node: 99 +Node: 100 +Node: 101 +Node: 102 +Node: 103 +Node: 104 +Node: 105 +Node: 106 +Node: 107 +Node: 108 +Node: 109 +Node: 110 +Node: 111 +Node: 112 +Node: 113 +Node: 114 +Node: 115 +Node: 116 +Node: 117 +Node: 118 +Node: 119 +Node: 120 +Node: 121 +Node: 122 +Node: 123 +Node: 124 +Node: 125 +Node: 126 +Node: 127 +Node: 128 +Node: 129 +Node: 130 +Node: 131 +Node: 132 +Node: 133 +Node: 134 +Node: 135 +Node: 136 +Node: 137 +Node: 138 +Node: 139 +Node: 140 +Node: 141 +Node: 142 +Node: 143 +Node: 144 +Node: 145 +Node: 146 +Node: 147 +Node: 148 +Node: 149 +Node: 150 +Node: 151 +Node: 152 +Node: 153 +Node: 154 +Node: 155 +Node: 156 +Node: 157 +Node: 158 +Node: 159 +Node: 160 +Node: 161 +Node: 162 +Node: 163 +Node: 164 +Node: 165 +Node: 166 +Node: 167 +Node: 168 +Node: 169 +Node: 170 +Node: 171 +Node: 172 +Node: 173 +Node: 174 +Node: 175 +Node: 176 +Node: 177 +Node: 178 +Node: 179 +Node: 180 +Node: 181 +Node: 182 +Node: 183 +Node: 184 +Node: 185 +Node: 186 +Node: 187 +Node: 188 +Node: 189 +Node: 190 +Node: 191 +Node: 192 +Node: 193 +Node: 194 +Node: 197 +Node: 198 +Node: 201 +Node: 202 +Node: 205 +Node: 206 +Node: 209 +Node: 210 +Node: 213 +Node: 214 +Node: 217 +Node: 218 +Node: 221 +Node: 222 +Node: 225 +Node: 226 +Node: 229 +Node: 230 +Node: 233 +Node: 234 +Node: 237 +Node: 238 + Group name: top +!#nodeset nodes +Shape. Dimension=0 +#Fields=0 +Node: 195 +Node: 196 +Node: 199 +Node: 200 +Node: 203 +Node: 204 +Node: 207 +Node: 208 +Node: 211 +Node: 212 +Node: 215 +Node: 216 +Node: 219 +Node: 220 +Node: 223 +Node: 224 +Node: 227 +Node: 228 +Node: 231 +Node: 232 +Node: 235 +Node: 236 +Node: 239 +Node: 240 +Node: 241 +Node: 242 +Node: 243 +Node: 244 +Node: 245 +Node: 246 +Node: 247 +Node: 248 +Node: 249 +Node: 250 +Node: 251 +Node: 252 +Node: 253 +Node: 254 +Node: 255 +Node: 256 +Node: 257 +Node: 258 +Node: 259 +Node: 260 +Node: 261 +Node: 262 +Node: 263 +Node: 264 +Node: 265 +Node: 266 +Node: 267 +Node: 268 +Node: 269 +Node: 270 +Node: 271 +Node: 272 +Node: 273 +Node: 274 +Node: 275 +Node: 276 +Node: 277 +Node: 278 +Node: 279 +Node: 280 +Node: 281 +Node: 282 +Node: 283 +Node: 284 +Node: 285 +Node: 286 +Node: 287 +Node: 288 diff --git a/tests/test_fitcube.py b/tests/test_fitcube.py new file mode 100644 index 0000000..56fe067 --- /dev/null +++ b/tests/test_fitcube.py @@ -0,0 +1,179 @@ +import math +import os +import unittest +from opencmiss.utils.zinc.field import createFieldMeshIntegral +from opencmiss.zinc.result import RESULT_OK +from scaffoldfitter.fitter import Fitter +from scaffoldfitter.fitterjson import decodeJSONFitterSteps +from scaffoldfitter.fitterstepalign import FitterStepAlign +from scaffoldfitter.fitterstepfit import FitterStepFit + +here = os.path.abspath(os.path.dirname(__file__)) + +def assertAlmostEqualList(testcase, actualList, expectedList, delta): + assert len(actualList) == len(expectedList) + for actual, expected in zip(actualList, expectedList): + testcase.assertAlmostEqual(actual, expected, delta=delta) + +def getRotationMatrix(eulerAngles): + """ + From OpenCMISS-Zinc graphics_library.cpp, transposed. + :param eulerAngles: 3-component field of angles in radians, components: + 1 = azimuth (about z) + 2 = elevation (about rotated y) + 3 = roll (about rotated x) + :return: 9-component rotation matrix varying fastest across, suitable for pre-multiplying [x, y, z]. + """ + cos_azimuth = math.cos(eulerAngles[0]) + sin_azimuth = math.sin(eulerAngles[0]) + cos_elevation = math.cos(eulerAngles[1]) + sin_elevation = math.sin(eulerAngles[1]) + cos_roll = math.cos(eulerAngles[2]) + sin_roll = math.sin(eulerAngles[2]) + return [ + cos_azimuth*cos_elevation, + cos_azimuth*sin_elevation*sin_roll - sin_azimuth*cos_roll, + cos_azimuth*sin_elevation*cos_roll + sin_azimuth*sin_roll, + sin_azimuth*cos_elevation, + sin_azimuth*sin_elevation*sin_roll + cos_azimuth*cos_roll, + sin_azimuth*sin_elevation*cos_roll - cos_azimuth*sin_roll, + -sin_elevation, + cos_elevation*sin_roll, + cos_elevation*cos_roll, + ] + +def transformCoordinatesList(xIn : list, transformationMatrix, translation): + """ + Transforms coordinates by multiplying by 9-component transformationMatrix + then offsetting by translation. + :xIn: List of 3-D coordinates to transform: + :return: List of 3-D transformed coordinates. + """ + assert (len(xIn) > 0) and (len(xIn[0]) == 3) and (len(transformationMatrix) == 9) and (len(translation) == 3) + xOut = [] + for x in xIn: + x2 = [] + for c in range(3): + v = translation[c] + for d in range(3): + v += transformationMatrix[c*3 + d]*x[d] + x2.append(v) + xOut.append(x2) + return xOut + +class FitCubeToSphereTestCase(unittest.TestCase): + + def test_alignFixedRandomData(self): + """ + Test alignment of model and data to known transformations. + """ + zinc_model_file = os.path.join(here, "resources", "cube_to_sphere.exf") + zinc_data_file = os.path.join(here, "resources", "cube_to_sphere_data_random.exf") + fitter = Fitter(zinc_model_file, zinc_data_file) + fitter.setDiagnosticLevel(1) + fitter.load() + + self.assertEqual(fitter.getModelCoordinatesField().getName(), "coordinates") + self.assertEqual(fitter.getDataCoordinatesField().getName(), "data_coordinates") + self.assertEqual(fitter.getMarkerGroup().getName(), "marker") + bottomCentre1 = fitter.evaluateNodeGroupMeanCoordinates("bottom", "coordinates", isData = False) + sidesCentre1 = fitter.evaluateNodeGroupMeanCoordinates("sides", "coordinates", isData = False) + topCentre1 = fitter.evaluateNodeGroupMeanCoordinates("top", "coordinates", isData = False) + assertAlmostEqualList(self, bottomCentre1, [ 0.5, 0.5, 0.0 ], delta=1.0E-7) + assertAlmostEqualList(self, sidesCentre1, [ 0.5, 0.5, 0.5 ], delta=1.0E-7) + assertAlmostEqualList(self, topCentre1, [ 0.5, 0.5, 1.0 ], delta=1.0E-7) + align = FitterStepAlign(fitter) + align.setScale(1.1) + align.setTranslation([ 0.1, -0.2, 0.3 ]) + align.setRotation([ math.pi/4.0, math.pi/8.0, math.pi/2.0 ]) + self.assertTrue(align.isAlignMarkers()) + align.setAlignMarkers(False) + align.run() + rotation = align.getRotation() + scale = align.getScale() + translation = align.getTranslation() + rotationMatrix = getRotationMatrix(rotation) + transformationMatrix = [ v*scale for v in rotationMatrix ] + bottomCentre2Expected, sidesCentre2Expected, topCentre2Expected = transformCoordinatesList( + [ bottomCentre1, sidesCentre1, topCentre1 ], transformationMatrix, translation) + bottomCentre2 = fitter.evaluateNodeGroupMeanCoordinates("bottom", "coordinates", isData = False) + sidesCentre2 = fitter.evaluateNodeGroupMeanCoordinates("sides", "coordinates", isData = False) + topCentre2 = fitter.evaluateNodeGroupMeanCoordinates("top", "coordinates", isData = False) + assertAlmostEqualList(self, bottomCentre2, bottomCentre2Expected, delta=1.0E-7) + assertAlmostEqualList(self, sidesCentre2, sidesCentre2Expected, delta=1.0E-7) + assertAlmostEqualList(self, topCentre2, topCentre2Expected, delta=1.0E-7) + + def test_alignMarkersFitRegularData(self): + """ + Test automatic alignment of model and data using fiducial markers. + """ + zinc_model_file = os.path.join(here, "resources", "cube_to_sphere.exf") + zinc_data_file = os.path.join(here, "resources", "cube_to_sphere_data_regular.exf") + fitter = Fitter(zinc_model_file, zinc_data_file) + fitter.setDiagnosticLevel(1) + fitter.load() + coordinates = fitter.getModelCoordinatesField() + self.assertEqual(coordinates.getName(), "coordinates") + self.assertEqual(fitter.getDataCoordinatesField().getName(), "data_coordinates") + self.assertEqual(fitter.getMarkerGroup().getName(), "marker") + #fitter.getRegion().writeFile(os.path.join(here, "resources", "km_fitgeometry1.exf")) + fieldmodule = fitter.getFieldmodule() + surfaceAreaField = createFieldMeshIntegral(coordinates, fitter.getMesh(2), number_of_points=4) + volumeField = createFieldMeshIntegral(coordinates, fitter.getMesh(3), number_of_points=3) + fieldcache = fieldmodule.createFieldcache() + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 6.0, delta=1.0E-6) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 1.0, delta=1.0E-7) + + align = FitterStepAlign(fitter) + self.assertTrue(align.isAlignMarkers()) + align.setAlignMarkers(True) + align.run() + #fitter.getRegion().writeFile(os.path.join(here, "resources", "km_fitgeometry2.exf")) + rotation = align.getRotation() + scale = align.getScale() + translation = align.getTranslation() + assertAlmostEqualList(self, rotation, [ -0.25*math.pi, 0.0, 0.0 ], delta=1.0E-4) + self.assertAlmostEqual(scale, 0.8047378476539072, places=5) + assertAlmostEqualList(self, translation, [ -0.5690355950594247, 1.1068454682130484e-05, -0.4023689233125251 ], delta=1.0E-6) + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 3.885618020657802, delta=1.0E-6) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 0.5211506471189844, delta=1.0E-6) + + fit1 = FitterStepFit(fitter) + fit1.setMarkerWeight(1.0) + fit1.setCurvaturePenaltyWeight(0.1) + fit1.setNumberOfIterations(3) + fit1.setUpdateReferenceState(True) + fit1.run() + #fitter.getRegion().writeFile(os.path.join(here, "resources", "km_fitgeometry3.exf")) + + result, surfaceArea = surfaceAreaField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(surfaceArea, 3.1892231780263853, delta=1.0E-4) + result, volume = volumeField.evaluateReal(fieldcache, 1) + self.assertEqual(result, RESULT_OK) + self.assertAlmostEqual(volume, 0.5276229458448985, delta=1.0E-4) + + # test json serialisation + s = fitter.encodeSettingsJSON() + fitter2 = Fitter(zinc_model_file, zinc_data_file) + fitter2.decodeSettingsJSON(s, decodeJSONFitterSteps) + fitterSteps = fitter2.getFitterSteps() + self.assertEqual(2, len(fitterSteps)) + self.assertTrue(isinstance(fitterSteps[0], FitterStepAlign)) + self.assertTrue(isinstance(fitterSteps[1], FitterStepFit)) + #fitter2.load() + #for fitterStep in fitterSteps: + # fitterStep.run() + s2 = fitter.encodeSettingsJSON() + self.assertEqual(s, s2) + +if __name__ == "__main__": + unittest.main()