diff --git a/README.rst b/README.rst index d7a96dd..f082a42 100644 --- a/README.rst +++ b/README.rst @@ -10,8 +10,10 @@ Building the example Instructions on how to configure and build with CMake:: git clone https://github.com/OpenCMISS-Examples/homogeneous_pipe_axial_extension.git + cd homogeneous_pipe_axial_extension mkdir build - cmake -DOpenCMISSLibs_DIR=/path/to/opencmisslib/install ../homogeneous_pipe_axial_extension + cd build + cmake -DOpenCMISS_INSTALL_ROOT=/path/to/opencmiss/install ../. make # cmake --build . will also work here and is much more platform agnostic. Running the example diff --git a/src/python/exfile.py b/src/python/exfile.py index 16fce57..b54eb12 100755 --- a/src/python/exfile.py +++ b/src/python/exfile.py @@ -1,662 +1,672 @@ -""" Module for reading exnode and exelem files used by cmgui - - Has only been used for reading exfiles produced by OpenCMISS - and doesn't attempt to be able to read any valid file. -""" - -import gzip -import numpy as np -import re - -__all__ = [ - "Exregion" - "Exnode", - "ExnodeField", - "ExnodeComponent", - "ExnodeNode", - "Exelem", - "ExelemField", - "ExelemComponent", - "ExelemElement", - "ExfileError", - ] - -class Exregion(object): - """ Store and retrieve data from an exelem file - """ - def __init__(self, filepath): - self.fields = [] - self.elements = [] - self.sections = [] - self.nodeids = [] - with FileWithLineNumber(filepath, 'r') as f: - self._read_header(f) - self.num_element_values = self._calc_num_element_values() - while True: - try: - self.sections.append(ExnodeSection(f, self)) - except ExfileError: - break -#Load the node ids - uid = set() - for section in self.sections: - for node in section.nodes: - uid.add(node.number) - self.nodeids = list(uid) -# Go back to the past 2 lines - f.rollbacktwice() -# read the exelem header data - self.num_dims = read_regex(f, r'Shape.\s+Dimension=\s*([0-9]+)') - self.num_scale_factor_sets = int(read_regex(f, r'#Scale factor sets=\s*([0-9]+)')) - self.num_scale_factors = 0 - for i in range(self.num_scale_factor_sets): - self.num_scale_factors += int(read_regex(f, r'#Scale factors\s*=\s*([0-9]+)')) - self.num_nodes = int(read_regex(f, r'#Nodes=\s*([0-9]+)')) - self.num_fields = int(read_regex(f, r'#Fields=\s*([0-9]+)')) - for i in range(self.num_fields): - field = ExelemField(f, self) - self.fields.append(field) - while True: - try: - self._read_element(f) - except EOFError: - break - self.num_elements = len(self.elements) - #self.num_nodes = sum(section.num_nodes for section in self.sections) - self.num_nodes = len(self.nodeids) - - def _read_header(self, f): - line = f.readline().strip() - regex = r'Group name|Region: ([A-Za-z0-9_\:\.\/]+)$' -#Skip root region specification - while True: - match = re.search(regex, line) - if match is not None: - if len(match.groups()) == 1: - self.group_name = match.group(1) - else: - self.group_name = match.groups() - line = f.readline().strip() - else: - break -#ZINC generated exregion files have !#nodeset nodes, check for that - regex = r'nodeset nodes$' - match = re.search(regex, line) - if match is None: - f.rollback(); - - def node_values(self, field_name, component_name, node_num): - """ Return all the field component derivative values - at the given node number - """ - for section in self.sections: - try: - return section.node_values(field_name, component_name, node_num) - except NodeNotFound: - pass - raise ValueError("Node %d not found in any exnode section." % node_num) - - def node_value(self, field_name, component_name, node_num, - derivative_number=1): - """ Return the field component value at the given node and derivative - Derivatives are numbered from 1, with 1 being no derivative. - """ - for section in self.sections: - try: - return section.node_value(field_name, component_name, node_num, - derivative_number) - except NodeNotFound: - pass - raise ValueError("Node %d not found in any exnode section." % node_num) - - - def _calc_num_element_values(self): - num_values = 0 - for field in self.fields: - for component in field.components: - if component.component_type == 'grid based': - component_num_values = np.product( - [i + 1 for i in component.divisions]) - num_values += component_num_values - return num_values - - def _read_element(self, f): - element_line = f.readline() - if element_line == "": - raise EOFError - try: - indices = map(int, element_line.split(':')[1].split()) - except: - print element_line - raise - if indices[1] == 0 and indices[2] == 0: - #raise ExfileError(f, "Face or line elements not supported") - values = [] - if self.num_element_values > 0: - expect_line(f, "Values:") - while len(values) < self.num_element_values: - line = f.readline() - values.extend(map(float, line.split())) - -#Ignore faces, lines - nodes = [] - element_line = f.readline() -#Move file pointer until Nodes are found - regex = r'Nodes:' - match = re.search(regex, element_line) - while match is None: - element_line = f.readline() - match = re.search(regex, element_line) - while len(nodes) < self.num_nodes: - line = f.readline() - nodes.extend(map(int, line.split())) - - scale_factors = [] - element_line = f.readline() - regex = r'Scale factors:' - if re.search(regex, element_line) is not None: - while len(scale_factors) < self.num_scale_factors: - line = f.readline() - scale_factors.extend(map(float, line.split())) - else: - f.rollback() - - self.elements.append(ExelemElement(indices, nodes, values, scale_factors)) - - def element_values(self, field_name, component_name, element_num): - """Return the all field component values at the given element number - """ - element = self.elements[element_num - 1] - value_index = 0 - for field in self.fields: - for component in field.components: - if component.component_type == 'grid based': - component_num_values = np.product( - [i + 1 for i in component.divisions]) - if (field.name == field_name and - component.name == str(component_name)): - return element.values[value_index:value_index + component_num_values] - else: - value_index += component_num_values - raise ValueError("Couldn't find field and component values") - - def quadPhi1(self,xi): - return 2*(xi-1.0)*(xi-0.5) - - def quadPhi2(self,xi): - return 4*xi*(1.0-xi) - - def quadPhi3(self,xi): - return 2*xi*(xi-0.5) - - def quad3Delembasisfactors(self,xi1,xi2,xi3): - functions = [self.quadPhi1,self.quadPhi2,self.quadPhi3] - phi=[] - for k in range(0,3): - phik = functions[k](xi3) - for j in range(0,3): - phij = functions[j](xi2) - for i in range(0,3): - phii = functions[i](xi1) - phi.append(phii*phij*phik) - return phi - -class Exnode(object): - """ Store and retrieve data from an exnode file - """ - def __init__(self, filepath): - self.sections = [] - with FileWithLineNumber(filepath, 'r') as f: - self._read_header(f) - while True: - try: - self.sections.append(ExnodeSection(f, self)) - except EOFError: - break - self.num_nodes = sum(section.num_nodes for section in self.sections) - - def _read_header(self, f): - self.group_name = read_regex(f, r'Group name|Region: ([A-Za-z0-9_\:\.]+)$') - - def node_values(self, field_name, component_name, node_num): - """ Return all the field component derivative values - at the given node number - """ - for section in self.sections: - try: - return section.node_values(field_name, component_name, node_num) - except NodeNotFound: - pass - raise ValueError("Node %d not found in any exnode section." % node_num) - - def node_value(self, field_name, component_name, node_num, - derivative_number=1): - """ Return the field component value at the given node and derivative - Derivatives are numbered from 1, with 1 being no derivative. - """ - for section in self.sections: - try: - return section.node_value(field_name, component_name, node_num, - derivative_number) - except NodeNotFound: - pass - raise ValueError("Node %d not found in any exnode section." % node_num) - - -class ExnodeSection(object): - def __init__(self, f, exnode): - self.exnode = exnode - self.fields = [] - self.nodes = [] - self._read_section_header(f) - self.num_node_values = self._calc_num_node_values() - while True: - try: - self._read_node(f) - except EOFError: - break - except ExfileError: - # Could be start of new section - f.rollback() - break - self.num_nodes = len(self.nodes) - - def _read_section_header(self, f): - try: - self.num_fields = int(read_regex(f, - r'#Fields=\s*([0-9]+)')) - except ExfileError: - if f.readline() == '': - raise EOFError - else: - raise - for i in range(self.num_fields): - field = ExnodeField(f, self) - self.fields.append(field) - - def _calc_num_node_values(self): - num_values = 0 - for field in self.fields: - for component in field.components: - num_values += 1 + component.num_derivatives - return num_values - - def _read_node(self, f): - line = f.readline().strip() - if line == "": - raise EOFError - number = int(read_string_regex(f, line, - r'Node:\s*([0-9]+)')) - read = 0 - values = np.empty(self.num_node_values) - while read < self.num_node_values: - line = f.readline() - try: - new_values = map(float, line.split()) - except ValueError: - raise ExfileError(f, "Expecting node values, got: %s" % line.strip()) - if read + len(new_values) > self.num_node_values: - raise ExfileError(f, "Got more node values than expected.") - values[read:read + len(new_values)] = new_values - read += len(new_values) - - self.nodes.append(ExnodeNode(number, values)) - - def _get_field_component(self, field_name, component_name): - try: - field = next(f for f in self.fields if f.name == field_name) - except StopIteration: - raise ValueError("Couldn't find field %s" % field_name) - try: - # Accept integer for component name - component_name = str(component_name) - component = next( - c for c in field.components if c.name == component_name) - except StopIteration: - raise ValueError("Couldn't find component %s" % component_name) - return field, component - - def node_values(self, field_name, component_name, node_num): - """ Return all the field component derivative values - at the given node number - """ - try: - node = next(n for n in self.nodes if n.number == node_num) - except StopIteration: - raise NodeNotFound() - field, component = self._get_field_component( - field_name, component_name) - - value_index = component.value_index - 1 - component_num_values = 1 + component.num_derivatives - return node.values[value_index:value_index + component_num_values] - - def node_value(self, field_name, component_name, node_num, derivative_number=1): - """ Return the field component value at the given node and derivative - Derivatives are numbered from 1, with 1 being no derivative. - """ - try: - node = next(n for n in self.nodes if n.number == node_num) - except StopIteration: - raise NodeNotFound() - field, component = self._get_field_component(field_name, component_name) - - value_index = component.value_index - 1 - value_index += derivative_number - 1 - if not 0 < derivative_number <= component.num_derivatives + 1: - raise ValueError("Invalid derivative number: %d" % derivative_number) - return node.values[value_index] - - -class ExnodeField(object): - """ A field definition from an exnode file - """ - def __init__(self, f, exnode): - self.exnode = exnode - - # Eg: - # 1) Geometry, coordinate, rectangular cartesian, #Components=3 - self.components = [] - - declaration = f.readline().split(',') - index, self.name = read_string_regex(f, declaration[0], - r'([0-9]+)\)\s+([A-Za-z0-9_\s\/]+)') - self.index = int(index) - self.num_components = int(read_string_regex(f, declaration[3], - r'#Components\s*=\s*([0-9]+)')) - for i in range(self.num_components): - component = ExnodeComponent(f, self) - self.components.append(component) - - def __repr__(self): - return '' % ( - self.index, self.name, self.num_components) - - -class ExnodeComponent(object): - """A field component definition from an exnode file - """ - def __init__(self, f, field): - self.field = field - # Eg: - # x. Value index= 1, #Derivatives= 0 - # or - # 1. Value index= 49, #Derivatives= 7(d/ds1,d/ds2,d2/ds1ds2,d/ds3,d2/ds1ds3,d2/ds2ds3,d3/ds1ds2ds3) - declaration = f.readline().strip().split(',', 1) - self.name = read_string_regex(f, declaration[0], - r'^([a-zA-Z0-9]+)\.') - self.value_index = int(read_string_regex(f, declaration[0], - r'Value index\s*=\s*([0-9]+)')) - self.num_derivatives = int(read_string_regex(f, declaration[1], - r'#Derivatives\s*=\s*([0-9]+)')) - if self.num_derivatives > 0: - derivative_name_list = read_string_regex(f, declaration[1], - r'#Derivatives\s*=\s*[0-9]+\s*\(([a-zA-Z0-9\,\/\s]+)\)') - self.derivative_names = [ - n.strip() for n in derivative_name_list.split(',')] - else: - self.derivative_names = [] - - def __repr__(self): - return '' % ( - self.name, self.num_derivatives) - - -class ExnodeNode(object): - """A node from an exnode file - """ - def __init__(self, number, values): - self.number = number - self.values = values - - def __repr__(self): - return "" % self.number - - -class Exelem(object): - """ Store and retrieve data from an exelem file - """ - def __init__(self, filepath): - self.fields = [] - self.elements = [] - with FileWithLineNumber(filepath, 'r') as f: - self._read_header(f) - self.num_element_values = self._calc_num_element_values() - while True: - try: - self._read_element(f) - except EOFError: - break - self.num_elements = len(self.elements) - - def _read_header(self, f): - self.group_name = read_regex(f, - r'Group name|Region: ([A-Za-z0-9_\:\.]+)') - self.num_dims = int(read_regex(f, - r'Shape.\s+Dimension=\s*([0-9]+)')) - self.num_scale_factor_sets = int(read_regex(f, - r'#Scale factor sets=\s*([0-9]+)')) - self.num_scale_factors = 0 - for i in range(self.num_scale_factor_sets): - self.num_scale_factors += int(read_regex(f, - r'#Scale factors\s*=\s*([0-9]+)')) - self.num_nodes = int(read_regex(f, - r'#Nodes=\s*([0-9]+)')) - self.num_fields = int(read_regex(f, - r'#Fields=\s*([0-9]+)')) - for i in range(self.num_fields): - field = ExelemField(f, self) - self.fields.append(field) - - def _calc_num_element_values(self): - num_values = 0 - for field in self.fields: - for component in field.components: - if component.component_type == 'grid based': - component_num_values = np.product( - [i + 1 for i in component.divisions]) - num_values += component_num_values - return num_values - - def _read_element(self, f): - element_line = f.readline() - if element_line == "": - raise EOFError - indices = map(int, element_line.split(':')[1].split()) - if indices[1] == 0 and indices[2] == 0: - # raise ExfileError(f, "Face or line elements not supported") - values = [] - if self.num_element_values > 0: - expect_line(f, "Values:") - while len(values) < self.num_element_values: - line = f.readline() - values.extend(map(float, line.split())) - - expect_line(f, "Nodes:") - nodes = [] - while len(nodes) < self.num_nodes: - line = f.readline() - nodes.extend(map(int, line.split())) - - expect_line(f, "Scale factors:") - scale_factors = [] - while len(scale_factors) < self.num_scale_factors: - line = f.readline() - scale_factors.extend(map(float, line.split())) - - self.elements.append( - ExelemElement(indices, nodes, values, scale_factors)) - - def element_values(self, field_name, component_name, element_num): - """Return the all field component values at the given element number - """ - element = self.elements[element_num - 1] - value_index = 0 - for field in self.fields: - for component in field.components: - if component.component_type == 'grid based': - component_num_values = np.product( - [i + 1 for i in component.divisions]) - if (field.name == field_name and - component.name == str(component_name)): - return element.values[value_index:value_index + component_num_values] - else: - value_index += component_num_values - raise ValueError("Couldn't find field and component values") - - -class ExelemField(object): - """A field definition from an exelem file - """ - def __init__(self, f, exelem): - self.exelem = exelem - - # Eg: - # 1) Geometry, coordinate, rectangular cartesian, #Components=3 - self.components = [] - - declaration = f.readline().split(',') - index, self.name = read_string_regex(f, declaration[0], - r'([0-9]+)\)\s+([A-Za-z0-9_\s\/]+)') - self.index = int(index) - self.num_components = int(read_string_regex(f, declaration[3], - r'#Components\s*=\s*([0-9]+)$')) - for i in range(self.num_components): - component = ExelemComponent(f, self) - self.components.append(component) - - def __repr__(self): - return '' % ( - self.index, self.name, self.num_components) - - -class ExelemComponent(object): - """A field component definition from an exelem file - """ - def __init__(self, f, field): - self.field = field - # Eg: - # x. l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based. - # #Nodes= 8 - declaration = f.readline().strip().split(',') - self.name = read_string_regex(f, declaration[0], r'^([a-zA-Z0-9]+)\.') - self.component_type = declaration[2].strip().strip('.') - if self.component_type == 'standard node based': - self._read_nodal_component(f) - elif self.component_type == 'grid based': - self._read_grid_component(f) - else: - raise ExfileError(f, "Unsupported component type: %s" % - self.component_type) - - def _read_nodal_component(self, f): - self.num_nodes = int(read_regex(f, r'#Nodes\s*=\s*([0-9]+)')) - self.node_num_values = {} - self.value_indices = {} - self.scale_factor_indices = {} - for node in range(1, self.num_nodes + 1): - self.node_num_values[node] = read_regex(f, - r'[0-9]+\.\s*#Values\s*=\s*([0-9]+)') - value_indices_line = f.readline().strip() - scale_factor_indices_line = f.readline().strip() - self.value_indices[node] = map(int, - value_indices_line.split(':')[1].strip().split()) - self.scale_factor_indices[node] = map(int, - scale_factor_indices_line.split(':')[1].strip().split()) - - def _read_grid_component(self, f): - grids = f.readline() - divisions = [d.strip().split('=')[1] for d in grids.split(',')] - self.divisions = map(int, divisions) - - def __repr__(self): - return '' % ( - self.component_type, self.name) - - -class ExelemElement(object): - """An element from an exelem file - """ - def __init__(self, indices, nodes, values, scale_factors): - self.indices = indices - self.number = indices[0] - self.nodes = nodes - self.values = values - self.scale_factors = scale_factors - - def __repr__(self): - return "" % self.indices[0] - - def __str__(self): - node_list = " ".join(str(n) for n in self.nodes) - return "%d: %s" % (self.indices[0], node_list) - - -class FileWithLineNumber(object): - def __init__ (self, path, *args): - if path.endswith('.gz'): - self.file = gzip.open(path, *args) - else: - self.file = open(path, *args) - self.linenum = 0 - self.prev_pos = self.file.tell() - self.tprev_pos, self.cur_pos = self.prev_pos, self.prev_pos - - def __enter__ (self): - return self - - def readline(self): - self.linenum += 1 - line = self.file.readline() - self.tprev_pos, self.prev_pos, self.cur_pos = self.prev_pos, self.cur_pos, self.file.tell() - return line - - def rollbacktwice(self): - self.file.seek(self.tprev_pos) - self.prev_pos = self.file.tell() - self.tprev_pos, self.cur_pos = self.prev_pos, self.prev_pos - - def rollback(self): - self.file.seek(self.prev_pos) - self.prev_pos = self.file.tell() - self.tprev_pos, self.cur_pos = self.prev_pos, self.prev_pos - - def __exit__ (self, exc_type, exc_value, traceback): - self.file.close() - - -class NodeNotFound(KeyError): - pass - - -class ExfileError(ValueError): - def __init__(self, file, message): - new_message = "Line %d: %s" % (file.linenum, message) - super(ExfileError, self).__init__(new_message) - - -def expect_line(f, expected): - line = f.readline().strip() - if line != expected: - raise ExfileError(f, "Expected '%s', got '%s'" % (expected, line)) - - -def read_regex(f, regex): - line = f.readline().strip() - match = re.search(regex, line) - if match is None: - raise ExfileError(f, "Expected '%s', got '%s'" % (regex, line)) - if len(match.groups()) == 1: - return match.group(1) - else: - return match.groups() - - -def read_string_regex(f, string, regex): - match = re.search(regex, string) - if match is None: - raise ExfileError(f, "Expected '%s', got '%s'" % (regex, string)) - if len(match.groups()) == 1: - return match.group(1) - else: - return match.groups() +""" Module for reading exnode and exelem files used by cmgui + + Has only been used for reading exfiles produced by OpenCMISS + and doesn't attempt to be able to read any valid file. +""" + +import gzip +import numpy as np +import re + +__all__ = [ + "Exregion" + "Exnode", + "ExnodeField", + "ExnodeComponent", + "ExnodeNode", + "Exelem", + "ExelemField", + "ExelemComponent", + "ExelemElement", + "ExfileError", + ] + +class Exregion(object): + """ Store and retrieve data from an exelem file + """ + def __init__(self, filepath): + self.fields = [] + self.elements = [] + self.sections = [] + self.nodeids = [] + with FileWithLineNumber(filepath, 'r') as f: + self._read_header(f) + self.num_element_values = self._calc_num_element_values() + while True: + try: + self.sections.append(ExnodeSection(f, self)) + except ExfileError: + break +#Load the node ids + uid = set() + for section in self.sections: + for node in section.nodes: + uid.add(node.number) + self.nodeids = list(uid) +# Go back to the past 2 lines + f.rollbacktwice() +# read the exelem header data + self.num_dims = read_regex(f, r'Shape.\s+Dimension=\s*([0-9]+)') + self.num_scale_factor_sets = int(read_regex(f, r'#Scale factor sets=\s*([0-9]+)')) + self.num_scale_factors = 0 + for i in range(self.num_scale_factor_sets): + self.num_scale_factors += int(read_regex(f, r'#Scale factors\s*=\s*([0-9]+)')) + self.num_nodes = int(read_regex(f, r'#Nodes=\s*([0-9]+)')) + self.num_fields = int(read_regex(f, r'#Fields=\s*([0-9]+)')) + for i in range(self.num_fields): + field = ExelemField(f, self) + self.fields.append(field) + while True: + try: + self._read_element(f) + except EOFError: + break + self.num_elements = len(self.elements) + #self.num_nodes = sum(section.num_nodes for section in self.sections) + self.num_nodes = len(self.nodeids) + + def _read_header(self, f): + line = f.readline().strip() + regex = r'Group name|Region: ([A-Za-z0-9_\:\.\/]+)$' +#Skip root region specification + while True: + match = re.search(regex, line) + if match is not None: + if len(match.groups()) == 1: + self.group_name = match.group(1) + else: + self.group_name = match.groups() + line = f.readline().strip() + else: + break +#ZINC generated exregion files have !#nodeset nodes, check for that + regex = r'nodeset nodes$' + match = re.search(regex, line) + if match is None: + f.rollback(); + + def node_values(self, field_name, component_name, node_num): + """ Return all the field component derivative values + at the given node number + """ + for section in self.sections: + try: + return section.node_values(field_name, component_name, node_num) + except NodeNotFound: + pass + raise ValueError("Node %d not found in any exnode section." % node_num) + + def node_value(self, field_name, component_name, node_num, + derivative_number=1): + """ Return the field component value at the given node and derivative + Derivatives are numbered from 1, with 1 being no derivative. + """ + for section in self.sections: + try: + return section.node_value(field_name, component_name, node_num, + derivative_number) + except NodeNotFound: + pass + raise ValueError("Node %d not found in any exnode section." % node_num) + + + def _calc_num_element_values(self): + num_values = 0 + for field in self.fields: + for component in field.components: + if component.component_type == 'grid based': + component_num_values = np.product( + [i + 1 for i in component.divisions]) + num_values += component_num_values + return num_values + + def _read_element(self, f): + element_line = f.readline() + if element_line == "": + raise EOFError + try: + indices = list(element_line.split(':')[1].split()) + except: + print(element_line) + raise + indices = str_list_to_int_list(indices) + if indices[1] == 0 and indices[2] == 0: + #raise ExfileError(f, "Face or line elements not supported") + values = [] + if self.num_element_values > 0: + expect_line(f, "Values:") + while len(values) < self.num_element_values: + line = f.readline() + values.extend(str_list_to_float_list(list(line.split()))) + +#Ignore faces, lines + nodes = [] + element_line = f.readline() +#Move file pointer until Nodes are found + regex = r'Nodes:' + match = re.search(regex, element_line) + while match is None: + element_line = f.readline() + match = re.search(regex, element_line) + while len(nodes) < self.num_nodes: + line = f.readline() + nodes.extend(str_list_to_int_list(list(line.split()))) + + scale_factors = [] + element_line = f.readline() + regex = r'Scale factors:' + if re.search(regex, element_line) is not None: + while len(scale_factors) < self.num_scale_factors: + line = f.readline() + scale_factors.extend(str_list_to_float_list(list(line.split()))) + else: + f.rollback() + + self.elements.append(ExelemElement(indices, nodes, values, scale_factors)) + + def element_values(self, field_name, component_name, element_num): + """Return the all field component values at the given element number + """ + element = self.elements[element_num - 1] + value_index = 0 + for field in self.fields: + for component in field.components: + if component.component_type == 'grid based': + component_num_values = np.product( + [i + 1 for i in component.divisions]) + if (field.name == field_name and + component.name == str(component_name)): + return element.values[value_index:value_index + component_num_values] + else: + value_index += component_num_values + raise ValueError("Couldn't find field and component values") + + def quadPhi1(self,xi): + return 2*(xi-1.0)*(xi-0.5) + + def quadPhi2(self,xi): + return 4*xi*(1.0-xi) + + def quadPhi3(self,xi): + return 2*xi*(xi-0.5) + + def quad3Delembasisfactors(self,xi1,xi2,xi3): + functions = [self.quadPhi1,self.quadPhi2,self.quadPhi3] + phi=[] + for k in range(0,3): + phik = functions[k](xi3) + for j in range(0,3): + phij = functions[j](xi2) + for i in range(0,3): + phii = functions[i](xi1) + phi.append(phii*phij*phik) + return phi + +class Exnode(object): + """ Store and retrieve data from an exnode file + """ + def __init__(self, filepath): + self.sections = [] + with FileWithLineNumber(filepath, 'r') as f: + self._read_header(f) + while True: + try: + self.sections.append(ExnodeSection(f, self)) + except EOFError: + break + self.num_nodes = sum(section.num_nodes for section in self.sections) + + def _read_header(self, f): + self.group_name = read_regex(f, r'Group name|Region: ([A-Za-z0-9_\:\.]+)$') + + def node_values(self, field_name, component_name, node_num): + """ Return all the field component derivative values + at the given node number + """ + for section in self.sections: + try: + return section.node_values(field_name, component_name, node_num) + except NodeNotFound: + pass + raise ValueError("Node %d not found in any exnode section." % node_num) + + def node_value(self, field_name, component_name, node_num, + derivative_number=1): + """ Return the field component value at the given node and derivative + Derivatives are numbered from 1, with 1 being no derivative. + """ + for section in self.sections: + try: + return section.node_value(field_name, component_name, node_num, + derivative_number) + except NodeNotFound: + pass + raise ValueError("Node %d not found in any exnode section." % node_num) + + +class ExnodeSection(object): + def __init__(self, f, exnode): + self.exnode = exnode + self.fields = [] + self.nodes = [] + self._read_section_header(f) + self.num_node_values = self._calc_num_node_values() + while True: + try: + self._read_node(f) + except EOFError: + break + except ExfileError: + # Could be start of new section + f.rollback() + break + self.num_nodes = len(self.nodes) + + def _read_section_header(self, f): + try: + self.num_fields = int(read_regex(f, + r'#Fields=\s*([0-9]+)')) + except ExfileError: + if f.readline() == '': + raise EOFError + else: + raise + for i in range(self.num_fields): + field = ExnodeField(f, self) + self.fields.append(field) + + def _calc_num_node_values(self): + num_values = 0 + for field in self.fields: + for component in field.components: + num_values += 1 + component.num_derivatives + return num_values + + def _read_node(self, f): + line = f.readline().strip() + if line == "": + raise EOFError + number = int(read_string_regex(f, line, + r'Node:\s*([0-9]+)')) + read = 0 + values = np.empty(self.num_node_values) + while read < self.num_node_values: + line = f.readline() + try: + new_values = list(line.split()) + except ValueError: + raise ExfileError(f, "Expecting node values, got: %s" % line.strip()) + str_list_to_float_list(new_values) + num_new_values = len(new_values) + if read + num_new_values > self.num_node_values: + raise ExfileError(f, "Got more node values than expected.") + values[read:read + num_new_values] = new_values + read += num_new_values + + self.nodes.append(ExnodeNode(number, values)) + + def _get_field_component(self, field_name, component_name): + try: + field = next(f for f in self.fields if f.name == field_name) + except StopIteration: + raise ValueError("Couldn't find field %s" % field_name) + try: + # Accept integer for component name + component_name = str(component_name) + component = next( + c for c in field.components if c.name == component_name) + except StopIteration: + raise ValueError("Couldn't find component %s" % component_name) + return field, component + + def node_values(self, field_name, component_name, node_num): + """ Return all the field component derivative values + at the given node number + """ + try: + node = next(n for n in self.nodes if n.number == node_num) + except StopIteration: + raise NodeNotFound() + field, component = self._get_field_component( + field_name, component_name) + + value_index = component.value_index - 1 + component_num_values = 1 + component.num_derivatives + return node.values[value_index:value_index + component_num_values] + + def node_value(self, field_name, component_name, node_num, derivative_number=1): + """ Return the field component value at the given node and derivative + Derivatives are numbered from 1, with 1 being no derivative. + """ + try: + node = next(n for n in self.nodes if n.number == node_num) + except StopIteration: + raise NodeNotFound() + field, component = self._get_field_component(field_name, component_name) + + value_index = component.value_index - 1 + value_index += derivative_number - 1 + if not 0 < derivative_number <= component.num_derivatives + 1: + raise ValueError("Invalid derivative number: %d" % derivative_number) + return node.values[value_index] + + +class ExnodeField(object): + """ A field definition from an exnode file + """ + def __init__(self, f, exnode): + self.exnode = exnode + + # Eg: + # 1) Geometry, coordinate, rectangular cartesian, #Components=3 + self.components = [] + + declaration = f.readline().split(',') + index, self.name = read_string_regex(f, declaration[0], + r'([0-9]+)\)\s+([A-Za-z0-9_\s\/]+)') + self.index = int(index) + self.num_components = int(read_string_regex(f, declaration[3], + r'#Components\s*=\s*([0-9]+)')) + for i in range(self.num_components): + component = ExnodeComponent(f, self) + self.components.append(component) + + def __repr__(self): + return '' % ( + self.index, self.name, self.num_components) + + +class ExnodeComponent(object): + """A field component definition from an exnode file + """ + def __init__(self, f, field): + self.field = field + # Eg: + # x. Value index= 1, #Derivatives= 0 + # or + # 1. Value index= 49, #Derivatives= 7(d/ds1,d/ds2,d2/ds1ds2,d/ds3,d2/ds1ds3,d2/ds2ds3,d3/ds1ds2ds3) + declaration = f.readline().strip().split(',', 1) + self.name = read_string_regex(f, declaration[0], + r'^([a-zA-Z0-9]+)\.') + self.value_index = int(read_string_regex(f, declaration[0], + r'Value index\s*=\s*([0-9]+)')) + self.num_derivatives = int(read_string_regex(f, declaration[1], + r'#Derivatives\s*=\s*([0-9]+)')) + if self.num_derivatives > 0: + derivative_name_list = read_string_regex(f, declaration[1], + r'#Derivatives\s*=\s*[0-9]+\s*\(([a-zA-Z0-9\,\/\s]+)\)') + self.derivative_names = [ + n.strip() for n in derivative_name_list.split(',')] + else: + self.derivative_names = [] + + def __repr__(self): + return '' % ( + self.name, self.num_derivatives) + + +class ExnodeNode(object): + """A node from an exnode file + """ + def __init__(self, number, values): + self.number = number + self.values = values + + def __repr__(self): + return "" % self.number + + +class Exelem(object): + """ Store and retrieve data from an exelem file + """ + def __init__(self, filepath): + self.fields = [] + self.elements = [] + with FileWithLineNumber(filepath, 'r') as f: + self._read_header(f) + self.num_element_values = self._calc_num_element_values() + while True: + try: + self._read_element(f) + except EOFError: + break + self.num_elements = len(self.elements) + + def _read_header(self, f): + self.group_name = read_regex(f, + r'Group name|Region: ([A-Za-z0-9_\:\.]+)') + self.num_dims = int(read_regex(f, + r'Shape.\s+Dimension=\s*([0-9]+)')) + self.num_scale_factor_sets = int(read_regex(f, + r'#Scale factor sets=\s*([0-9]+)')) + self.num_scale_factors = 0 + for i in range(self.num_scale_factor_sets): + self.num_scale_factors += int(read_regex(f, + r'#Scale factors\s*=\s*([0-9]+)')) + self.num_nodes = int(read_regex(f, + r'#Nodes=\s*([0-9]+)')) + self.num_fields = int(read_regex(f, + r'#Fields=\s*([0-9]+)')) + for i in range(self.num_fields): + field = ExelemField(f, self) + self.fields.append(field) + + def _calc_num_element_values(self): + num_values = 0 + for field in self.fields: + for component in field.components: + if component.component_type == 'grid based': + component_num_values = np.product( + [i + 1 for i in component.divisions]) + num_values += component_num_values + return num_values + + def _read_element(self, f): + element_line = f.readline() + if element_line == "": + raise EOFError + indices = list(element_line.split(':')[1].split()) + indices = str_list_to_int_list(indices) + if indices[1] == 0 and indices[2] == 0: + # raise ExfileError(f, "Face or line elements not supported") + values = [] + if self.num_element_values > 0: + expect_line(f, "Values:") + while len(values) < self.num_element_values: + line = f.readline() + values.extend(str_list_to_float_list(list(line.split()))) + + expect_line(f, "Nodes:") + nodes = [] + while len(nodes) < self.num_nodes: + line = f.readline() + nodes.extend(str_list_to_int_list(list(line.split()))) + + expect_line(f, "Scale factors:") + scale_factors = [] + while len(scale_factors) < self.num_scale_factors: + line = f.readline() + scale_factors.extend(str_list_to_float_list(list(line.split()))) + + self.elements.append( + ExelemElement(indices, nodes, values, scale_factors)) + + def element_values(self, field_name, component_name, element_num): + """Return the all field component values at the given element number + """ + element = self.elements[element_num - 1] + value_index = 0 + for field in self.fields: + for component in field.components: + if component.component_type == 'grid based': + component_num_values = np.product( + [i + 1 for i in component.divisions]) + if (field.name == field_name and + component.name == str(component_name)): + return element.values[value_index:value_index + component_num_values] + else: + value_index += component_num_values + raise ValueError("Couldn't find field and component values") + + +class ExelemField(object): + """A field definition from an exelem file + """ + def __init__(self, f, exelem): + self.exelem = exelem + + # Eg: + # 1) Geometry, coordinate, rectangular cartesian, #Components=3 + self.components = [] + + declaration = f.readline().split(',') + index, self.name = read_string_regex(f, declaration[0], + r'([0-9]+)\)\s+([A-Za-z0-9_\s\/]+)') + self.index = int(index) + self.num_components = int(read_string_regex(f, declaration[3], + r'#Components\s*=\s*([0-9]+)$')) + for i in range(self.num_components): + component = ExelemComponent(f, self) + self.components.append(component) + + def __repr__(self): + return '' % ( + self.index, self.name, self.num_components) + + +class ExelemComponent(object): + """A field component definition from an exelem file + """ + def __init__(self, f, field): + self.field = field + # Eg: + # x. l.Lagrange*l.Lagrange*l.Lagrange, no modify, standard node based. + # #Nodes= 8 + declaration = f.readline().strip().split(',') + self.name = read_string_regex(f, declaration[0], r'^([a-zA-Z0-9]+)\.') + self.component_type = declaration[2].strip().strip('.') + if self.component_type == 'standard node based': + self._read_nodal_component(f) + elif self.component_type == 'grid based': + self._read_grid_component(f) + else: + raise ExfileError(f, "Unsupported component type: %s" % + self.component_type) + + def _read_nodal_component(self, f): + self.num_nodes = int(read_regex(f, r'#Nodes\s*=\s*([0-9]+)')) + self.node_num_values = {} + self.value_indices = {} + self.scale_factor_indices = {} + for node in range(1, self.num_nodes + 1): + self.node_num_values[node] = read_regex(f, + r'[0-9]+\.\s*#Values\s*=\s*([0-9]+)') + value_indices_line = f.readline().strip() + scale_factor_indices_line = f.readline().strip() + self.value_indices[node] = str_list_to_int_list(list(value_indices_line.split(':')[1].strip().split())) + self.scale_factor_indices[node] = str_list_to_int_list(list(scale_factor_indices_line.split(':')[1].strip().split())) + + def _read_grid_component(self, f): + grids = f.readline() + divisions = [d.strip().split('=')[1] for d in grids.split(',')] + self.divisions = str_list_to_int_list(list(divisions)) + + def __repr__(self): + return '' % ( + self.component_type, self.name) + + +class ExelemElement(object): + """An element from an exelem file + """ + def __init__(self, indices, nodes, values, scale_factors): + self.indices = indices + self.number = indices[0] + self.nodes = nodes + self.values = values + self.scale_factors = scale_factors + + def __repr__(self): + return "" % self.indices[0] + + def __str__(self): + node_list = " ".join(str(n) for n in self.nodes) + return "%d: %s" % (self.indices[0], node_list) + + +class FileWithLineNumber(object): + def __init__ (self, path, *args): + if path.endswith('.gz'): + self.file = gzip.open(path, *args) + else: + self.file = open(path, *args) + self.linenum = 0 + self.prev_pos = self.file.tell() + self.tprev_pos, self.cur_pos = self.prev_pos, self.prev_pos + + def __enter__ (self): + return self + + def readline(self): + self.linenum += 1 + line = self.file.readline() + self.tprev_pos, self.prev_pos, self.cur_pos = self.prev_pos, self.cur_pos, self.file.tell() + return line + + def rollbacktwice(self): + self.file.seek(self.tprev_pos) + self.prev_pos = self.file.tell() + self.tprev_pos, self.cur_pos = self.prev_pos, self.prev_pos + + def rollback(self): + self.file.seek(self.prev_pos) + self.prev_pos = self.file.tell() + self.tprev_pos, self.cur_pos = self.prev_pos, self.prev_pos + + def __exit__ (self, exc_type, exc_value, traceback): + self.file.close() + + +class NodeNotFound(KeyError): + pass + + +class ExfileError(ValueError): + def __init__(self, file, message): + new_message = "Line %d: %s" % (file.linenum, message) + super(ExfileError, self).__init__(new_message) + + +def expect_line(f, expected): + line = f.readline().strip() + if line != expected: + raise ExfileError(f, "Expected '%s', got '%s'" % (expected, line)) + + +def read_regex(f, regex): + line = f.readline().strip() + match = re.search(regex, line) + if match is None: + raise ExfileError(f, "Expected '%s', got '%s'" % (regex, line)) + if len(match.groups()) == 1: + return match.group(1) + else: + return match.groups() + + +def read_string_regex(f, string, regex): + match = re.search(regex, string) + if match is None: + raise ExfileError(f, "Expected '%s', got '%s'" % (regex, string)) + if len(match.groups()) == 1: + return match.group(1) + else: + return match.groups() + +def str_list_to_int_list(str_list): + int_list = [int(n) for n in str_list] + return int_list + +def str_list_to_float_list(str_list): + float_list = [float(n) for n in str_list] + return float_list diff --git a/src/python/homogeneous_pipe_axial_extension.py b/src/python/homogeneous_pipe_axial_extension.py index a0e4d0c..c7b37f7 100755 --- a/src/python/homogeneous_pipe_axial_extension.py +++ b/src/python/homogeneous_pipe_axial_extension.py @@ -16,7 +16,7 @@ numberOfXi = 3 # Intialise OpenCMISS -from opencmiss.iron import iron +from opencmiss.opencmiss import OpenCMISS_Python as oc # Set problem parameters #Use pressure to enforce incompressibililty constraint @@ -25,6 +25,7 @@ numberOfGaussXi = 3 #Setup field number handles +contextUserNumber = 1 coordinateSystemUserNumber = 1 regionUserNumber = 1 basisUserNumber = 1 @@ -32,6 +33,7 @@ meshUserNumber = 1 cellMLUserNumber = 1 decompositionUserNumber = 1 +decomposerUserNumber = 1 equationsSetUserNumber = 1 problemUserNumber = 1 #Mesh component numbers @@ -46,76 +48,89 @@ cellMLIntermediateFieldUserNumber = 6 equationsSetFieldUserNumber = 7 +#quit() + # Set all diganostic levels on for testing -#iron.DiagnosticsSetOn(iron.DiagnosticTypes.ALL,[1,2,3,4,5],"Diagnostics",["DOMAIN_MAPPINGS_LOCAL_FROM_GLOBAL_CALCULATE"]) +#oc.DiagnosticsSetOn(oc.DiagnosticTypes.ALL,[1,2,3,4,5],"Diagnostics",["DOMAIN_MAPPINGS_LOCAL_FROM_GLOBAL_CALCULATE"]) if(usePressureBasis): numberOfMeshComponents = 2 else: numberOfMeshComponents = 1 +context = oc.Context() +context.Create(contextUserNumber) + +worldRegion = oc.Region() +context.WorldRegionGet(worldRegion) + # Set the OpenCMISS random seed so that we can test this example by using the # same parallel decomposition -numberOfRandomSeeds = iron.RandomSeedsSizeGet() +numberOfRandomSeeds = context.RandomSeedsSizeGet() randomSeeds = [0]*numberOfRandomSeeds randomSeeds[0] = 100 -iron.RandomSeedsSet(randomSeeds) +context.RandomSeedsSet(randomSeeds) # Get the number of computational nodes and this computational node number -numberOfComputationalNodes = iron.ComputationalNumberOfNodesGet() -computationalNodeNumber = iron.ComputationalNodeNumberGet() +computationEnvironment = oc.ComputationEnvironment() +context.ComputationEnvironmentGet(computationEnvironment) + +worldWorkGroup = oc.WorkGroup() +computationEnvironment.WorldWorkGroupGet(worldWorkGroup) +numberOfComputationalNodes = worldWorkGroup.NumberOfGroupNodesGet() +computationalNodeNumber = worldWorkGroup.GroupNodeNumberGet() if computationalNodeNumber == 0: if not os.path.exists("./results"): os.makedirs("./results") # Create a 3D rectangular cartesian coordinate system -coordinateSystem = iron.CoordinateSystem() -coordinateSystem.CreateStart(coordinateSystemUserNumber) +coordinateSystem = oc.CoordinateSystem() +coordinateSystem.CreateStart(coordinateSystemUserNumber,context) coordinateSystem.DimensionSet(3) coordinateSystem.CreateFinish() # Create a region and assign the coordinate system to the region -region = iron.Region() -region.CreateStart(regionUserNumber,iron.WorldRegion) +region = oc.Region() +region.CreateStart(regionUserNumber,worldRegion) region.LabelSet("Region") region.coordinateSystem = coordinateSystem region.CreateFinish() # Define basis -basis = iron.Basis() -basis.CreateStart(basisUserNumber) -basis.TypeSet(iron.BasisTypes.LAGRANGE_HERMITE_TP) +basis = oc.Basis() +basis.CreateStart(basisUserNumber,context) +basis.TypeSet(oc.BasisTypes.LAGRANGE_HERMITE_TP) basis.NumberOfXiSet(numberOfXi) -basis.InterpolationXiSet([iron.BasisInterpolationSpecifications.QUADRATIC_LAGRANGE]*numberOfXi) +basis.InterpolationXiSet([oc.BasisInterpolationSpecifications.QUADRATIC_LAGRANGE]*numberOfXi) if(numberOfGaussXi>0): basis.QuadratureNumberOfGaussXiSet([numberOfGaussXi]*numberOfXi) basis.CreateFinish() if(usePressureBasis): # Define pressure basis - pressureBasis = iron.Basis() - pressureBasis.CreateStart(pressureBasisUserNumber) - pressureBasis.TypeSet(iron.BasisTypes.LAGRANGE_HERMITE_TP) + pressureBasis = oc.Basis() + pressureBasis.CreateStart(pressureBasisUserNumber,context) + pressureBasis.TypeSet(oc.BasisTypes.LAGRANGE_HERMITE_TP) pressureBasis.NumberOfXiSet(numberOfXi) - pressureBasis.InterpolationXiSet([iron.BasisInterpolationSpecifications.LINEAR_LAGRANGE]*numberOfXi) + pressureBasis.InterpolationXiSet([oc.BasisInterpolationSpecifications.LINEAR_LAGRANGE]*numberOfXi) if(numberOfGaussXi>0): pressureBasis.QuadratureNumberOfGaussXiSet([numberOfGaussXi]*numberOfXi) pressureBasis.CreateFinish() # Start the creation of input mesh in the region -mesh = iron.Mesh() +mesh = oc.Mesh() mesh.CreateStart(meshUserNumber, region, numberOfXi) mesh.NumberOfComponentsSet(numberOfMeshComponents) mesh.NumberOfElementsSet(exregion.num_elements) # Define nodes for the mesh -nodes = iron.Nodes() +nodes = oc.Nodes() nodes.CreateStart(region, exregion.num_nodes) nodes.CreateFinish() #Specify the elementwise topology -quadraticElemnts = iron.MeshElements() +quadraticElemnts = oc.MeshElements() quadraticElemnts.CreateStart(mesh, quadraticMeshComponentNumber, basis) for elem in exregion.elements: quadraticElemnts.NodesSet(elem.number, elem.nodes) @@ -123,7 +138,7 @@ if(usePressureBasis): #Specify the elementwise topology that suits the mesh element basis - linearElements = iron.MeshElements() + linearElements = oc.MeshElements() linearElements.CreateStart(mesh, linearMeshComponentNumber, pressureBasis) for elem in exregion.elements: linearNodes = list(elem.nodes[i - 1] for i in [1, 3, 7, 9, 19, 21, 25, 27]) @@ -133,21 +148,25 @@ # Create a decomposition for the mesh -decomposition = iron.Decomposition() +decomposition = oc.Decomposition() decomposition.CreateStart(decompositionUserNumber,mesh) -decomposition.TypeSet(iron.DecompositionTypes.CALCULATED) -decomposition.NumberOfDomainsSet(numberOfComputationalNodes) decomposition.CreateFinish() +# Decompose +decomposer = oc.Decomposer() +decomposer.CreateStart(decomposerUserNumber,worldRegion,worldWorkGroup) +decompositionIndex = decomposer.DecompositionAdd(decomposition) +decomposer.CreateFinish() + # Create a field for the geometry -geometricField = iron.Field() +geometricField = oc.Field() geometricField.CreateStart(geometricFieldUserNumber,region) -geometricField.MeshDecompositionSet(decomposition) -geometricField.TypeSet(iron.FieldTypes.GEOMETRIC) -geometricField.VariableLabelSet(iron.FieldVariableTypes.U,"coordinates") -geometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,1,1) -geometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,2,1) -geometricField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,3,1) +geometricField.DecompositionSet(decomposition) +geometricField.TypeSet(oc.FieldTypes.GEOMETRIC) +geometricField.VariableLabelSet(oc.FieldVariableTypes.U,"coordinates") +geometricField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,1,1) +geometricField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,2,1) +geometricField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,3,1) geometricField.CreateFinish() #Set the geometric information from the exregion file @@ -157,7 +176,7 @@ # DOC-START define node coordinates # Read the geometric field for nodeNumber in exregion.nodeids: - nodeDomain = decomposition.NodeDomainGet(nodeNumber,1) + nodeDomain = decomposition.NodeDomainGet(1,nodeNumber) if(nodeDomain == computationalNodeNumber): version = 1 derivative = 1 @@ -166,8 +185,8 @@ component_name = ["1", "2", "3"][component - 1] value = exregion.node_value("coordinates", component_name, nodeNumber, derivative) geometricField.ParameterSetUpdateNodeDP( - iron.FieldVariableTypes.U, - iron.FieldParameterSetTypes.VALUES, + oc.FieldVariableTypes.U, + oc.FieldParameterSetTypes.VALUES, version, derivative, nodeNumber, component, value) coord.append(value) #The cylinder has an axial length of 5 units and coord[2] contains the axial coordinate value @@ -178,113 +197,113 @@ # Update the geometric field parameters -geometricField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES) -geometricField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES) +geometricField.ParameterSetUpdateStart(oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES) +geometricField.ParameterSetUpdateFinish(oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES) # DOC-END define node coordinates # Create a fibre field and attach it to the geometric field -fibreField = iron.Field() +fibreField = oc.Field() fibreField.CreateStart(fibreFieldUserNumber,region) -fibreField.TypeSet(iron.FieldTypes.FIBRE) -fibreField.MeshDecompositionSet(decomposition) +fibreField.TypeSet(oc.FieldTypes.FIBRE) +fibreField.DecompositionSet(decomposition) fibreField.GeometricFieldSet(geometricField) -fibreField.VariableLabelSet(iron.FieldVariableTypes.U,"Fibre") +fibreField.VariableLabelSet(oc.FieldVariableTypes.U,"Fibre") fibreField.CreateFinish() # Create the dependent field #Create the dependent field with 4 variables and the respective number of components #1 U_Var_Type 4 components: 3 displacement (quad interpol) + 1 pressure (lin interpol)) -#2 DELUDELN_Var_Type 4 components: 3 displacement (quad interpol) + 1 pressure (lin interpol)) +#2 T_Var_Type 4 components: 3 displacement (quad interpol) + 1 pressure (lin interpol)) #3 U1_Var_Type 6 components: 6 independent components of the strain tensor (quad interpol) [independent] #4 U2_Var_Type 6 components: 6 independent components of the stress tensor (quad interpol) [dependent] -dependentField = iron.Field() +dependentField = oc.Field() dependentField.CreateStart(dependentFieldUserNumber,region) -dependentField.VariableLabelSet(iron.FieldVariableTypes.U,"Dependent") -dependentField.TypeSet(iron.FieldTypes.GEOMETRIC_GENERAL) -dependentField.MeshDecompositionSet(decomposition) +dependentField.VariableLabelSet(oc.FieldVariableTypes.U,"Dependent") +dependentField.TypeSet(oc.FieldTypes.GEOMETRIC_GENERAL) +dependentField.DecompositionSet(decomposition) dependentField.GeometricFieldSet(geometricField) -dependentField.DependentTypeSet(iron.FieldDependentTypes.DEPENDENT) +dependentField.DependentTypeSet(oc.FieldDependentTypes.DEPENDENT) #Displacement, Gradient, Strain and Stress Tensors dependentField.NumberOfVariablesSet(4) -dependentField.VariableTypesSet([iron.FieldVariableTypes.U,iron.FieldVariableTypes.DELUDELN,iron.FieldVariableTypes.U1,iron.FieldVariableTypes.U2]) -dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U,4) -dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.DELUDELN,4) -dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U1,6) -dependentField.NumberOfComponentsSet(iron.FieldVariableTypes.U2,6) +dependentField.VariableTypesSet([oc.FieldVariableTypes.U,oc.FieldVariableTypes.T,oc.FieldVariableTypes.U1,oc.FieldVariableTypes.U2]) +dependentField.NumberOfComponentsSet(oc.FieldVariableTypes.U,4) +dependentField.NumberOfComponentsSet(oc.FieldVariableTypes.T,4) +dependentField.NumberOfComponentsSet(oc.FieldVariableTypes.U1,6) +dependentField.NumberOfComponentsSet(oc.FieldVariableTypes.U2,6) #Assign the mesh from which the quantities are determined -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,1,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,2,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,3,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,1,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,2,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,3,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,1,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,2,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,3,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T,1,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T,2,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T,3,quadraticMeshComponentNumber) if(usePressureBasis): - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,4,linearMeshComponentNumber) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,4,linearMeshComponentNumber) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,4,linearMeshComponentNumber) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T,4,linearMeshComponentNumber) else: - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,4,quadraticMeshComponentNumber) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,4,quadraticMeshComponentNumber) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,4,quadraticMeshComponentNumber) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T,4,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1,1,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1,2,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1,3,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1,4,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1,5,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U1,6,quadraticMeshComponentNumber) - -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2,1,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2,2,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2,3,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2,4,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2,5,quadraticMeshComponentNumber) -dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U2,6,quadraticMeshComponentNumber) - -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,1,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,2,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,3,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,4,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,5,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U1,6,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) - -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,1,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,2,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,3,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,4,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,5,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) -dependentField.ComponentInterpolationSet(iron.FieldVariableTypes.U2,6,iron.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1,1,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1,2,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1,3,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1,4,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1,5,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U1,6,quadraticMeshComponentNumber) + +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2,1,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2,2,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2,3,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2,4,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2,5,quadraticMeshComponentNumber) +dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U2,6,quadraticMeshComponentNumber) + +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U1,1,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U1,2,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U1,3,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U1,4,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U1,5,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U1,6,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) + +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U2,1,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U2,2,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U2,3,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U2,4,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U2,5,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) +dependentField.ComponentInterpolationSet(oc.FieldVariableTypes.U2,6,oc.FieldInterpolationTypes.GAUSS_POINT_BASED) if(usePressureBasis): - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.U,4,2) - dependentField.ComponentMeshComponentSet(iron.FieldVariableTypes.DELUDELN,4,2) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.U,4,2) + dependentField.ComponentMeshComponentSet(oc.FieldVariableTypes.T,4,2) -dependentField.VariableLabelSet(iron.FieldVariableTypes.U1,"strain") -dependentField.VariableLabelSet(iron.FieldVariableTypes.U2,"stress") +dependentField.VariableLabelSet(oc.FieldVariableTypes.U1,"strain") +dependentField.VariableLabelSet(oc.FieldVariableTypes.U2,"stress") dependentField.CreateFinish() # Initialise dependent field from undeformed geometry and displacement bcs and set hydrostatic pressure -iron.Field.ParametersToFieldParametersComponentCopy( - geometricField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1, - dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,1) -iron.Field.ParametersToFieldParametersComponentCopy( - geometricField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,2, - dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,2) -iron.Field.ParametersToFieldParametersComponentCopy( - geometricField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,3, - dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,3) +oc.Field.ParametersToFieldParametersComponentCopy( + geometricField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,1, + dependentField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,1) +oc.Field.ParametersToFieldParametersComponentCopy( + geometricField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,2, + dependentField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,2) +oc.Field.ParametersToFieldParametersComponentCopy( + geometricField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,3, + dependentField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,3) #Set hydrostatic pressure, the value depends on the constitutive law being used -iron.Field.ComponentValuesInitialiseDP(dependentField,iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES,4,-8.0) +oc.Field.ComponentValuesInitialiseDP(dependentField,oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES,4,-8.0) # Create the equations_set -equationsSetField = iron.Field() -equationsSet = iron.EquationsSet() -equationsSetSpecification = [iron.EquationsSetClasses.ELASTICITY, - iron.EquationsSetTypes.FINITE_ELASTICITY, - iron.EquationsSetSubtypes.CONSTITUTIVE_LAW_IN_CELLML_EVALUATE] +equationsSetField = oc.Field() +equationsSet = oc.EquationsSet() +equationsSetSpecification = [oc.EquationsSetClasses.ELASTICITY, + oc.EquationsSetTypes.FINITE_ELASTICITY, + oc.EquationsSetSubtypes.CONSTITUTIVE_LAW_IN_CELLML_EVALUATE] equationsSet.CreateStart(equationsSetUserNumber, region, fibreField, equationsSetSpecification, equationsSetFieldUserNumber, equationsSetField) equationsSet.CreateFinish() @@ -295,7 +314,7 @@ #DOC-START create CellML environment # Create the CellML environment -cellML = iron.CellML() +cellML = oc.CellML() cellML.CreateStart(cellMLUserNumber, region) # Import a Mooney-Rivlin material law from a file mooneyRivlinModel = cellML.ModelImport("mooney_rivlin.xml") @@ -303,19 +322,19 @@ #DOC-START flag variables # Now we have imported the model we are able to specify which variables from the model we want to set from openCMISS -cellML.VariableSetAsKnown(mooneyRivlinModel, "equations/E11") -cellML.VariableSetAsKnown(mooneyRivlinModel, "equations/E12") -cellML.VariableSetAsKnown(mooneyRivlinModel, "equations/E13") -cellML.VariableSetAsKnown(mooneyRivlinModel, "equations/E22") -cellML.VariableSetAsKnown(mooneyRivlinModel, "equations/E23") -cellML.VariableSetAsKnown(mooneyRivlinModel, "equations/E33") +cellML.VariableSetAsKnown(mooneyRivlinModel, "main/E11") +cellML.VariableSetAsKnown(mooneyRivlinModel, "main/E12") +cellML.VariableSetAsKnown(mooneyRivlinModel, "main/E13") +cellML.VariableSetAsKnown(mooneyRivlinModel, "main/E22") +cellML.VariableSetAsKnown(mooneyRivlinModel, "main/E23") +cellML.VariableSetAsKnown(mooneyRivlinModel, "main/E33") # and variables to get from the CellML -cellML.VariableSetAsWanted(mooneyRivlinModel, "equations/Tdev11") -cellML.VariableSetAsWanted(mooneyRivlinModel, "equations/Tdev12") -cellML.VariableSetAsWanted(mooneyRivlinModel, "equations/Tdev13") -cellML.VariableSetAsWanted(mooneyRivlinModel, "equations/Tdev22") -cellML.VariableSetAsWanted(mooneyRivlinModel, "equations/Tdev23") -cellML.VariableSetAsWanted(mooneyRivlinModel, "equations/Tdev33") +cellML.VariableSetAsWanted(mooneyRivlinModel, "main/Tdev11") +cellML.VariableSetAsWanted(mooneyRivlinModel, "main/Tdev12") +cellML.VariableSetAsWanted(mooneyRivlinModel, "main/Tdev13") +cellML.VariableSetAsWanted(mooneyRivlinModel, "main/Tdev22") +cellML.VariableSetAsWanted(mooneyRivlinModel, "main/Tdev23") +cellML.VariableSetAsWanted(mooneyRivlinModel, "main/Tdev33") #DOC-END flag variables #DOC-START create cellml finish @@ -328,22 +347,22 @@ #Now we can set up the field variable component <--> CellML model variable mappings. #DOC-START map strain components #Map the strain components -cellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,1, iron.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"equations/E11", iron.FieldParameterSetTypes.VALUES) -cellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,2, iron.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"equations/E12", iron.FieldParameterSetTypes.VALUES) -cellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,3, iron.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"equations/E13", iron.FieldParameterSetTypes.VALUES) -cellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,4, iron.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"equations/E22", iron.FieldParameterSetTypes.VALUES) -cellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,5, iron.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"equations/E23", iron.FieldParameterSetTypes.VALUES) -cellML.CreateFieldToCellMLMap(dependentField,iron.FieldVariableTypes.U1,6, iron.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"equations/E33", iron.FieldParameterSetTypes.VALUES) +cellML.CreateFieldToCellMLMap(dependentField,oc.FieldVariableTypes.U1,1, oc.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"main/E11", oc.FieldParameterSetTypes.VALUES) +cellML.CreateFieldToCellMLMap(dependentField,oc.FieldVariableTypes.U1,2, oc.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"main/E12", oc.FieldParameterSetTypes.VALUES) +cellML.CreateFieldToCellMLMap(dependentField,oc.FieldVariableTypes.U1,3, oc.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"main/E13", oc.FieldParameterSetTypes.VALUES) +cellML.CreateFieldToCellMLMap(dependentField,oc.FieldVariableTypes.U1,4, oc.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"main/E22", oc.FieldParameterSetTypes.VALUES) +cellML.CreateFieldToCellMLMap(dependentField,oc.FieldVariableTypes.U1,5, oc.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"main/E23", oc.FieldParameterSetTypes.VALUES) +cellML.CreateFieldToCellMLMap(dependentField,oc.FieldVariableTypes.U1,6, oc.FieldParameterSetTypes.VALUES,mooneyRivlinModel,"main/E33", oc.FieldParameterSetTypes.VALUES) #DOC-END map strain components #DOC-START map stress components #Map the stress components -cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"equations/Tdev11", iron.FieldParameterSetTypes.VALUES,dependentField,iron.FieldVariableTypes.U2,1,iron.FieldParameterSetTypes.VALUES) -cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"equations/Tdev12", iron.FieldParameterSetTypes.VALUES,dependentField,iron.FieldVariableTypes.U2,2,iron.FieldParameterSetTypes.VALUES) -cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"equations/Tdev13", iron.FieldParameterSetTypes.VALUES,dependentField,iron.FieldVariableTypes.U2,3,iron.FieldParameterSetTypes.VALUES) -cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"equations/Tdev22", iron.FieldParameterSetTypes.VALUES,dependentField,iron.FieldVariableTypes.U2,4,iron.FieldParameterSetTypes.VALUES) -cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"equations/Tdev23", iron.FieldParameterSetTypes.VALUES,dependentField,iron.FieldVariableTypes.U2,5,iron.FieldParameterSetTypes.VALUES) -cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"equations/Tdev33", iron.FieldParameterSetTypes.VALUES,dependentField,iron.FieldVariableTypes.U2,6,iron.FieldParameterSetTypes.VALUES) +cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"main/Tdev11", oc.FieldParameterSetTypes.VALUES,dependentField,oc.FieldVariableTypes.U2,1,oc.FieldParameterSetTypes.VALUES) +cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"main/Tdev12", oc.FieldParameterSetTypes.VALUES,dependentField,oc.FieldVariableTypes.U2,2,oc.FieldParameterSetTypes.VALUES) +cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"main/Tdev13", oc.FieldParameterSetTypes.VALUES,dependentField,oc.FieldVariableTypes.U2,3,oc.FieldParameterSetTypes.VALUES) +cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"main/Tdev22", oc.FieldParameterSetTypes.VALUES,dependentField,oc.FieldVariableTypes.U2,4,oc.FieldParameterSetTypes.VALUES) +cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"main/Tdev23", oc.FieldParameterSetTypes.VALUES,dependentField,oc.FieldVariableTypes.U2,5,oc.FieldParameterSetTypes.VALUES) +cellML.CreateCellMLToFieldMap(mooneyRivlinModel,"main/Tdev33", oc.FieldParameterSetTypes.VALUES,dependentField,oc.FieldVariableTypes.U2,6,oc.FieldParameterSetTypes.VALUES) #DOC-END map stress components #Finish the creation of cellML <--> OpenCMISS field maps @@ -351,7 +370,7 @@ #DOC-START define CellML models field #Create the CellML models field -cellMLModelsField = iron.Field() +cellMLModelsField = oc.Field() cellML.ModelsFieldCreateStart(cellMLModelsFieldUserNumber,cellMLModelsField) cellML.ModelsFieldCreateFinish() @@ -372,108 +391,108 @@ for xk in range(0,numberOfGaussXi): xi3 = (1.0+xk)*xidiv gaussPointIdx = gaussPointIdx + 1 - cellMLModelsField.ParameterSetUpdateGaussPoint(iron.FieldVariableTypes.U, - iron.FieldParameterSetTypes.VALUES, + cellMLModelsField.ParameterSetUpdateGaussPoint(oc.FieldVariableTypes.U, + oc.FieldParameterSetTypes.VALUES, gaussPointIdx,elment.number,1, mooneyRivlinModel) # Update the models field parameters -cellMLModelsField.ParameterSetUpdateStart(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES) -cellMLModelsField.ParameterSetUpdateFinish(iron.FieldVariableTypes.U,iron.FieldParameterSetTypes.VALUES) +cellMLModelsField.ParameterSetUpdateStart(oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES) +cellMLModelsField.ParameterSetUpdateFinish(oc.FieldVariableTypes.U,oc.FieldParameterSetTypes.VALUES) #DOC-END define CellML models field #DOC-START define CellML parameters and intermediate fields #Create the CellML parameters field --- the strain field -cellMLParametersField = iron.Field() +cellMLParametersField = oc.Field() cellML.ParametersFieldCreateStart(cellMLParametersFieldUserNumber,cellMLParametersField) cellML.ParametersFieldCreateFinish() # Create the CellML intermediate field --- the stress field -cellMLIntermediateField = iron.Field() +cellMLIntermediateField = oc.Field() cellML.IntermediateFieldCreateStart(cellMLIntermediateFieldUserNumber,cellMLIntermediateField) cellML.IntermediateFieldCreateFinish() #DOC-END define CellML parameters and intermediate fields # Create equations -equations = iron.Equations() +equations = oc.Equations() equationsSet.EquationsCreateStart(equations) -equations.SparsityTypeSet(iron.EquationsSparsityTypes.SPARSE) -equations.OutputTypeSet(iron.EquationsOutputTypes.NONE) +equations.SparsityTypeSet(oc.EquationsSparsityTypes.SPARSE) +equations.OutputTypeSet(oc.EquationsOutputTypes.NONE) equationsSet.EquationsCreateFinish() #DOC-START define CellML finite elasticity problem #Define the problem -problem = iron.Problem() -problemSpecification = [iron.ProblemClasses.ELASTICITY, - iron.ProblemTypes.FINITE_ELASTICITY, - iron.ProblemSubtypes.FINITE_ELASTICITY_WITH_CELLML] -problem.CreateStart(problemUserNumber, problemSpecification) +problem = oc.Problem() +problemSpecification = [oc.ProblemClasses.ELASTICITY, + oc.ProblemTypes.FINITE_ELASTICITY, + oc.ProblemSubtypes.FINITE_ELASTICITY_WITH_CELLML] +problem.CreateStart(problemUserNumber,context,problemSpecification) problem.CreateFinish() #DOC-END define CellML finite elasticity problem #Create the problem control loop problem.ControlLoopCreateStart() -ControlLoop = iron.ControlLoop() -problem.ControlLoopGet([iron.ControlLoopIdentifiers.NODE],ControlLoop) -ControlLoop.TypeSet(iron.ProblemControlLoopTypes.SIMPLE) +ControlLoop = oc.ControlLoop() +problem.ControlLoopGet([oc.ControlLoopIdentifiers.NODE],ControlLoop) +ControlLoop.TypeSet(oc.ControlLoopTypes.SIMPLE) problem.ControlLoopCreateFinish() #Create the problem solvers -nonLinearSolver = iron.Solver() -linearSolver = iron.Solver() +nonLinearSolver = oc.Solver() +linearSolver = oc.Solver() problem.SolversCreateStart() -problem.SolverGet([iron.ControlLoopIdentifiers.NODE],1,nonLinearSolver) -nonLinearSolver.OutputTypeSet(iron.SolverOutputTypes.MONITOR) -nonLinearSolver.NewtonJacobianCalculationTypeSet(iron.JacobianCalculationTypes.FD) +problem.SolverGet([oc.ControlLoopIdentifiers.NODE],1,nonLinearSolver) +nonLinearSolver.OutputTypeSet(oc.SolverOutputTypes.MONITOR) +nonLinearSolver.NewtonJacobianCalculationTypeSet(oc.JacobianCalculationTypes.EQUATIONS) nonLinearSolver.NewtonLinearSolverGet(linearSolver) #Use the DIRECT MUMPS solver -linearSolver.LinearTypeSet(iron.LinearSolverTypes.DIRECT) +linearSolver.LinearTypeSet(oc.LinearSolverTypes.DIRECT) #For large problems or problems with complex material behaviour, the direct solver may fail #In such cases either preconditioners or the following solvers can be tried #In case the matrix has zeros for some rows, SUPERLU is a good solver to try as it will report such errors -#linearSolver.LibraryTypeSet(iron.SolverLibraries.PASTIX) -#linearSolver.LibraryTypeSet(iron.SolverLibraries.SUPERLU) +#linearSolver.LibraryTypeSet(oc.SolverLibraries.PASTIX) +#linearSolver.LibraryTypeSet(oc.SolverLibraries.SUPERLU) problem.SolversCreateFinish() #DOC-START define CellML solver #Create the problem solver CellML equations -cellMLSolver = iron.Solver() +cellMLSolver = oc.Solver() problem.CellMLEquationsCreateStart() nonLinearSolver.NewtonCellMLSolverGet(cellMLSolver) -cellMLEquations = iron.CellMLEquations() +cellMLEquations = oc.CellMLEquations() cellMLSolver.CellMLEquationsGet(cellMLEquations) cellMLEquations.CellMLAdd(cellML) problem.CellMLEquationsCreateFinish() #DOC-END define CellML solver #Create the problem solver equations -solver = iron.Solver() -solverEquations = iron.SolverEquations() +solver = oc.Solver() +solverEquations = oc.SolverEquations() problem.SolverEquationsCreateStart() -problem.SolverGet([iron.ControlLoopIdentifiers.NODE],1,solver) +problem.SolverGet([oc.ControlLoopIdentifiers.NODE],1,solver) solver.SolverEquationsGet(solverEquations) -solverEquations.SparsityTypeSet(iron.SolverEquationsSparsityTypes.SPARSE) +solverEquations.SparsityTypeSet(oc.SolverEquationsSparsityTypes.SPARSE) equationsSetIndex = solverEquations.EquationsSetAdd(equationsSet) problem.SolverEquationsCreateFinish() # Prescribe boundary conditions (absolute nodal parameters) -boundaryConditions = iron.BoundaryConditions() +boundaryConditions = oc.BoundaryConditions() solverEquations.BoundaryConditionsCreateStart(boundaryConditions) #Here we model axial stretch, by pulling along the axis at both the ends of the cylinder and holding them (Direchlet) for nodeNumber in left_boundary_nodes: - nodeDomain = decomposition.NodeDomainGet(nodeNumber, 1) + nodeDomain = decomposition.NodeDomainGet(1,nodeNumber) if nodeDomain == computationalNodeNumber : - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, 1, nodeNumber, 1, iron.BoundaryConditionsTypes.FIXED, 0.0) - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, 1, nodeNumber, 2, iron.BoundaryConditionsTypes.FIXED, 0.0) - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, 1, nodeNumber, 3, iron.BoundaryConditionsTypes.FIXED, -1.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, 1, nodeNumber, 1, oc.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, 1, nodeNumber, 2, oc.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, 1, nodeNumber, 3, oc.BoundaryConditionsTypes.FIXED, -1.0) for nodeNumber in right_boundary_nodes: - nodeDomain = decomposition.NodeDomainGet(nodeNumber, 1) + nodeDomain = decomposition.NodeDomainGet(1,nodeNumber) if nodeDomain == computationalNodeNumber : - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, 1, nodeNumber, 1, iron.BoundaryConditionsTypes.FIXED, 0.0) - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, 1, nodeNumber, 2, iron.BoundaryConditionsTypes.FIXED, 0.0) - boundaryConditions.AddNode(dependentField, iron.FieldVariableTypes.U, 1, 1, nodeNumber, 3, iron.BoundaryConditionsTypes.FIXED, 1.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, 1, nodeNumber, 1, oc.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, 1, nodeNumber, 2, oc.BoundaryConditionsTypes.FIXED, 0.0) + boundaryConditions.AddNode(dependentField, oc.FieldVariableTypes.U, 1, 1, nodeNumber, 3, oc.BoundaryConditionsTypes.FIXED, 1.0) solverEquations.BoundaryConditionsCreateFinish() @@ -487,7 +506,7 @@ problem.Solve() # Export the results, here we export them as standard exnode, exelem files -fields = iron.Fields() +fields = oc.Fields() fields.CreateRegion(region) fields.NodesExport("./results/AxialStretch","FORTRAN") fields.ElementsExport("./results/AxialStretch","FORTRAN") diff --git a/src/python/mooney_rivlin.xml b/src/python/mooney_rivlin.xml index 9e02ec2..7723e2f 100755 --- a/src/python/mooney_rivlin.xml +++ b/src/python/mooney_rivlin.xml @@ -1,311 +1,154 @@ - - - - - - - - - - - Nickerson - David - - - d.nickerson@auckland.ac.nz - - - - The University of Auckland - The Bioengineering Institute - - - - - 2003-11-28 - - - - - - - - This is a CellML version of the Mooney-Rivlin constitutive material law, - defining the relation between the eight independent strain components - and the stress components. It is assumed that the strain components - will be controlled externally by the application using this CellML - model. - - - - pubmed_id - - - - - - - Master - Andre - T - - - - - - - Bob - Billy - - - - - - - What cool article to reference ?? - - - year - - - The Journal of Cool Stuff - - volume - 1 - 1000 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - We'll use this component as the "interface" to the model, all - other components are hidden via encapsulation in this component. - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - In this simple model we only have one component, which holds the - six equations. - - - - - - - - - - - - - - - - - - - - - - - - - - - - Tdev11 - - - 2.0 - c1 - - - 4.0 - c2 - - E22 - E33 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Tdev11 + + + + + 2.0 + c1 + + + + 4.0 + c2 + + + E22 + E33 + + + + + 4.0 + c2 + + - - - 4 - c2 - - - - - - - Tdev22 - - - 2.0 - c1 - - - 4.0 - c2 - - E11 - E33 + + + Tdev22 + + + + + 2.0 + c1 + + + + 4.0 + c2 + + + E11 + E33 + + + + + 4.0 + c2 + + - - - 4 - c2 - - - - - - - Tdev33 - - - 2.0 - c1 - - - 4.0 - c2 - - E11 - E22 + + + Tdev33 + + + + + 2.0 + c1 + + + + 4.0 + c2 + + + E11 + E22 + + + + + 4.0 + c2 + + - - - 4 - c2 - - - - - - - Tdev12 - - - 4.0 - E12 - c2 - - - - - - - Tdev13 - - - 4.0 - E13 - c2 - - - - - - - Tdev23 - - - 4.0 - E23 - c2 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + Tdev12 + + + + + 4.0 + + E12 + c2 + + + + + Tdev13 + + + + + 4.0 + + E13 + c2 + + + + + Tdev23 + + + + + 4.0 + + E23 + c2 + + + + +