From 26a8c1b6e37f8afeeca250c122ee771f86dfacdd Mon Sep 17 00:00:00 2001 From: joergbuchwald Date: Thu, 2 Feb 2023 01:29:55 +0100 Subject: [PATCH 1/6] add parameters dictionary --- ogs6py/classes/parameter_type.py | 138 +++++++++++++++++++++++++++++++ ogs6py/classes/parameters.py | 103 ++++++++++++++++++++++- ogs6py/ogs.py | 59 +++++++------ 3 files changed, 273 insertions(+), 27 deletions(-) create mode 100644 ogs6py/classes/parameter_type.py diff --git a/ogs6py/classes/parameter_type.py b/ogs6py/classes/parameter_type.py new file mode 100644 index 0000000..42f6084 --- /dev/null +++ b/ogs6py/classes/parameter_type.py @@ -0,0 +1,138 @@ +# -*- coding: utf-8 -*- +""" +Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +from fastcore.utils import * +from lxml import etree as ET + + +# pylint: disable=C0103, R0902, R0914, R0913 + +class Parameter_type: + def __init__(self, xmlobject:object=None) -> None: + self.__dict__ = {} + self.xmlobject = xmlobject + if not self.xmlobject is None: + for parameter_property in self.xmlobject: + if parameter_property.tag == "expression": + try: + self.__dict__["expression"].append(parameter_property.text) + except: + self.__dict__["expression"] = [] + self.__dict__["expression"].append(parameter_property.text) + else: + self.__dict__[parameter_property.tag] = parameter_property.text + def __setitem__(self, key, item): + if not key in self.__dict__: + raise RuntimeError("property is not existing") + if key == "type": + raise RuntimeError("The Type can't be changed.") + expression_counter = -1 + for parameter_property in self.xmlobject: + if parameter_property.tag == key: + if key == "expression": + expression_counter += 1 + parameter_property.text = str(item[expression_counter]) + else: + parameter_property.text = (item) + + def __getitem__(self, key): + if not (key == "xmlobject"): + return self.__dict__[key] + + def __repr__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return repr(newdict) + + def __len__(self): + return len(self.__dict__) + + def __delitem__(self, key): + pass + #del self.__dict__[key] + + def clear(self): + return self.__dict__.clear() + + def copy(self): + return self.__dict__.copy() + + def has_key(self, k): + if not (k == "xmlobject"): + return k in self.__dict__ + + def update(self, *args, **kwargs): + pass + # return self.__dict__.update(*args, **kwargs) + + def keys(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return newdict.keys() + + def values(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return newdict.values() + + def items(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return newdict.items() + + def pop(self, *args): + pass + #return self.__dict__.pop(*args) + + def __cmp__(self, dict_): + return self.__cmp__(self.__dict__, dict_) + + def __contains__(self, item): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return item in newdict + + def __iter__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return iter(newdict) + +class Constant(Parameter_type): + pass + +class Function(Parameter_type): + pass + +class MeshNode(Parameter_type): + pass + +class MeshElement(Parameter_type): + pass + +class CurveScaled(Parameter_type): + pass + +class TimeDependentHeterogeneousParameter(Parameter_type): + pass + +class RandomFieldMeshElementParameter(Parameter_type): + pass +class Group(Parameter_type): + pass \ No newline at end of file diff --git a/ogs6py/classes/parameters.py b/ogs6py/classes/parameters.py index be3e139..d9664eb 100644 --- a/ogs6py/classes/parameters.py +++ b/ogs6py/classes/parameters.py @@ -8,12 +8,13 @@ """ # pylint: disable=C0103, R0902, R0914, R0913 from ogs6py.classes import build_tree +from ogs6py.classes import parameter_type class Parameters(build_tree.BuildTree): """ Class for managing the parameters section of the project file. """ - def __init__(self): + def __init__(self, xmlobject=None): self.tree = { 'parameters': { 'tag': 'parameters', @@ -22,6 +23,106 @@ def __init__(self): 'children': {} } } + self.parameter = {} + self.xmlobject = xmlobject + #self.__dict__ = {} + if not (xmlobject is None): + for prmt in xmlobject: + for parameter_property in prmt: + if parameter_property.tag == "type": + param_type = parameter_property.text + elif parameter_property.tag == "name": + param_name = parameter_property.text + if param_type == "Constant": + self.__dict__[param_name] = parameter_type.Constant(prmt) + elif param_type == "Function": + self.__dict__[param_name] = parameter_type.Function(prmt) + elif param_type == "MeshNode": + self.__dict__[param_name] = parameter_type.MeshNode(prmt) + elif param_type == "MeshElement": + self.__dict__[param_name] = parameter_type.MeshElement(prmt) + elif param_type == "CurveScaled": + self.__dict__[param_name] = parameter_type.CurveScaled(prmt) + elif param_type == "TimeDependentHeterogeneousParameter": + self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) + elif param_type == "RandomFieldMeshElementParameter": + self.__dict__[param_name] = parameter_type.RandomFieldMeshElementParameter(prmt) + elif param_type == "Group": + self.__dict__[param_name] = parameter_type.Group(prmt) + + def __getitem__(self, key): + if not (key in ["tree","parameter", "xmlobject"]): + return self.__dict__[key] + + def __repr__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","parameter", "xmlobject"]): + newdict[k] = v + return repr(newdict) + + def __len__(self): + return len(self.__dict__) + + def __delitem__(self, key): + pass + #del self.__dict__[key] + + def clear(self): + return self.__dict__.clear() + + def copy(self): + return self.__dict__.copy() + + def has_key(self, k): + if not (k in ["tree","parameter", "xmlobject"]): + return k in self.__dict__ + + def update(self, *args, **kwargs): + pass + # return self.__dict__.update(*args, **kwargs) + + def keys(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","parameter", "xmlobject"]): + newdict[k] = v + return newdict.keys() + + def values(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","parameter", "xmlobject"]): + newdict[k] = v + return newdict.values() + + def items(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","parameter", "xmlobject"]): + newdict[k] = v + return newdict.items() + + def pop(self, *args): + pass + #return self.__dict__.pop(*args) + + def __cmp__(self, dict_): + return self.__cmp__(self.__dict__, dict_) + + def __contains__(self, item): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","parameter", "xmlobject"]): + newdict[k] = v + return item in newdict + + def __iter__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","parameter", "xmlobject"]): + newdict[k] = v + return iter(newdict) def add_parameter(self, **args): """ diff --git a/ogs6py/ogs.py b/ogs6py/ogs.py index 5aca83b..90d5306 100644 --- a/ogs6py/ogs.py +++ b/ogs6py/ogs.py @@ -18,6 +18,7 @@ import shutil import pandas as pd from lxml import etree as ET +from fastcore.utils import * from ogs6py.classes import (display, geo, mesh, python_script, processes, media, timeloop, local_coordinate_system, parameters, curves, processvars, linsolvers, nonlinsolvers) import ogs6py.log_parser.log_parser as parser @@ -26,22 +27,22 @@ class OGS: """Class for an OGS6 model. - In this class everything for an OGS5 model can be specified. + In this class everything for an OGS6 model can be specified. Parameters ---------- - PROJECT_FILE : `str`, optional + project_file : `str`, optional Filename of the output project file Default: default.prj - INPUT_FILE : `str`, optional + input_file : `str`, optional Filename of the input project file - XMLSTRING : `str`,optional - OMP_NUM_THREADS : `int`, optional + xmlstring : `str`,optional + omp_num_threads : `int`, optional Sets the environmentvariable before OGS execution to restrict number of OMP Threads - VERBOSE : `bool`, optional + verbose : `bool`, optional Default: False """ - def __init__(self, **args): + def __init__(self, input_file=None, project_file="default.prj", xmlstring=None, verbose=False, omp_num_threads=None, **args): self.geo = geo.Geo() self.mesh = mesh.Mesh() self.pyscript = python_script.PythonScript() @@ -49,7 +50,7 @@ def __init__(self, **args): self.media = media.Media() self.timeloop = timeloop.TimeLoop() self.local_coordinate_system = local_coordinate_system.LocalCoordinateSystem() - self.parameters = parameters.Parameters() + #self.parameters = parameters.Parameters() self.curves = curves.Curves() self.processvars = processvars.ProcessVars() self.linsolvers = linsolvers.LinSolvers() @@ -61,32 +62,38 @@ def __init__(self, **args): self.include_elements = [] self.include_files = [] self.add_includes = [] + store_attr() + # **args only fror backwards compatibility if "VERBOSE" in args: self.verbose = args["VERBOSE"] else: self.verbose = False if "OMP_NUM_THREADS" in args: - self.threads = args["OMP_NUM_THREADS"] + self.omp_num_threads = args["OMP_NUM_THREADS"] else: - self.threads = None + self.omp_num_threads = None if "PROJECT_FILE" in args: - self.prjfile = args['PROJECT_FILE'] + self.project_file = args['PROJECT_FILE'] else: print("PROJECT_FILE for output not given. Calling it default.prj.") - self.prjfile = "default.prj" + self.project_file = "default.prj" if "INPUT_FILE" in args: - if os.path.isfile(args['INPUT_FILE']) is True: - self.inputfile = args['INPUT_FILE'] + self.input_file = args['INPUT_FILE'] + if self.input_file is not None: + if os.path.isfile(self.input_file) is True: _ = self._get_root() if self.verbose is True: display.Display(self.tree) else: raise RuntimeError(f"Input project file {args['INPUT_FILE']} not found.") - else: - self.inputfile = None - if "XMLSTRING" in args: + if "XMLSTRING" in args or (self.xmlstring is not None): root = ET.fromstring(args['XMLSTRING']) self.tree = ET.ElementTree(root) + try: + paramobj = self.tree.find("./parameters") + except AttributeError: + paramobj = None + self.parameters = parameters.Parameters(xmlobject=paramobj) def __dict2xml(self, parent, dictionary): for entry in dictionary: @@ -115,8 +122,8 @@ def __replace_blocks_by_includes(self): def _get_root(self): if self.tree is None: - if self.inputfile is not None: - self.tree = ET.parse(self.inputfile) + if self.input_file is not None: + self.tree = ET.parse(self.input_file) else: self.build_tree() root = self.tree.getroot() @@ -566,10 +573,10 @@ def run_model(self, logfile="out.log", path=None, args=None, container_path=None """ ogs_path = "" - if self.threads is None: + if self.omp_num_threads is None: env_export = "" else: - env_export = f"export OMP_NUM_THREADS={self.threads} && " + env_export = f"export OMP_NUM_THREADS={self.omp_num_threads} && " if not container_path is None: container_path = os.path.expanduser(container_path) if os.path.isfile(container_path) is False: @@ -610,9 +617,9 @@ def run_model(self, logfile="out.log", path=None, args=None, container_path=None if not args is None: cmd += f"{args} " if write_logs is True: - cmd += f"{self.prjfile} > {self.logfile}" + cmd += f"{self.project_file} > {self.logfile}" else: - cmd += f"{self.prjfile}" + cmd += f"{self.project_file}" startt = time.time() if sys.platform == "win32": returncode = subprocess.run(cmd, shell=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT) @@ -621,7 +628,7 @@ def run_model(self, logfile="out.log", path=None, args=None, container_path=None stopt = time.time() self.exec_time = stopt - startt if returncode.returncode == 0: - print(f"OGS finished with project file {self.prjfile}.") + print(f"OGS finished with project file {self.project_file}.") print(f"Execution took {self.exec_time} s") else: print(f"Error code: {returncode.returncode}") @@ -679,7 +686,7 @@ def write_input(self, keep_includes=False): ET.indent(self.tree, space=" ") if self.verbose is True: display.Display(self.tree) - self.tree.write(self.prjfile, + self.tree.write(self.project_file, encoding="ISO-8859-1", xml_declaration=True, pretty_print=True) @@ -688,7 +695,7 @@ def write_input(self, keep_includes=False): ET.indent(self.tree, space=" ") if self.verbose is True: display.Display(self.tree) - self.tree.write(self.prjfile, + self.tree.write(self.project_file, encoding="ISO-8859-1", xml_declaration=True, pretty_print=True) From 3db0b4db98853c9a9f395e90bf5288cd9219b4bd Mon Sep 17 00:00:00 2001 From: joergbuchwald Date: Fri, 3 Feb 2023 01:17:08 +0100 Subject: [PATCH 2/6] add curve dict --- ogs6py/classes/curve.py | 111 ++++++++++++++++++++++++++++ ogs6py/classes/curves.py | 121 ++++++++++++++++++++++++++++++- ogs6py/classes/parameter_type.py | 97 ++++++++++++++++++++++--- ogs6py/classes/parameters.py | 83 +++++++++++++++++++-- ogs6py/ogs.py | 39 +++++++--- setup.py | 1 + 6 files changed, 421 insertions(+), 31 deletions(-) create mode 100644 ogs6py/classes/curve.py diff --git a/ogs6py/classes/curve.py b/ogs6py/classes/curve.py new file mode 100644 index 0000000..418ce11 --- /dev/null +++ b/ogs6py/classes/curve.py @@ -0,0 +1,111 @@ +# -*- coding: utf-8 -*- +""" +Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) + Distributed under a Modified BSD License. + See accompanying file LICENSE or + http://www.opengeosys.org/project/license + +""" +from fastcore.utils import * +import numpy as np +from lxml import etree as ET +# pylint: disable=C0103, R0902, R0914, R0913 + +class Curve: + def __init__(self, xmlobject:object=None) -> None: + self.__dict__ = {} + self.xmlobject = xmlobject + if not self.xmlobject is None: + for curve_property in self.xmlobject: + if curve_property.tag == "name": + self.__dict__[curve_property.tag] = curve_property.text + else: + self.__dict__[curve_property.tag] = curve_property.text.split(" ") + assert(len(self.__dict__["coords"])==len(self.__dict__["values"])) + + def __setitem__(self, key, item): + if not key in self.__dict__: + raise RuntimeError("property is not existing") + expression_counter = -1 + if key in ["coords", "values"]: + for i, entry in enumerate(item): + item[i] = str(entry) + for curve_property in self.xmlobject: + if curve_property.tag == key: + curve_property.text = ' '.join(item) + else: + for curve_property in self.xmlobject: + if curve_property.tag == key: + curve_property.text = item + + def __getitem__(self, key): + if not (key == "xmlobject"): + return self.__dict__[key] + + def __repr__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return repr(newdict) + + def __len__(self): + return len(self.__dict__) + + def __delitem__(self, key): + pass + #del self.__dict__[key] + + def clear(self): + return self.__dict__.clear() + + def copy(self): + return self.__dict__.copy() + + def has_key(self, k): + if not (k == "xmlobject"): + return k in self.__dict__ + + def update(self, *args, **kwargs): + pass + # return self.__dict__.update(*args, **kwargs) + + def keys(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return newdict.keys() + + def items(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return newdict.items() + + def pop(self, *args): + pass + #return self.__dict__.pop(*args) + + def __cmp__(self, dict_): + return self.__cmp__(self.__dict__, dict_) + + def __contains__(self, item): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return item in newdict + + def __iter__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k == "xmlobject"): + newdict[k] = v + return iter(newdict) + + def evaluate_values(self,x): + coords = np.array([float(val) for val in self.__dict__["coords"]]) + values = np.array([float(val) for val in self.__dict__["values"]]) + return np.interp(x, coords, values) diff --git a/ogs6py/classes/curves.py b/ogs6py/classes/curves.py index 287d8de..567d948 100644 --- a/ogs6py/classes/curves.py +++ b/ogs6py/classes/curves.py @@ -1,19 +1,21 @@ # -*- coding: utf-8 -*- """ -Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org) +Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org) Distributed under a Modified BSD License. See accompanying file LICENSE or http://www.opengeosys.org/project/license """ # pylint: disable=C0103, R0902, R0914, R0913 +from lxml import etree as ET from ogs6py.classes import build_tree +from ogs6py.classes import curve class Curves(build_tree.BuildTree): """ Class to create the curve section of the project file. """ - def __init__(self): + def __init__(self, xmlobject=None): self.tree = { 'curves': { 'tag': 'curves', @@ -22,6 +24,121 @@ def __init__(self): 'children': {} } } + self.xmlobject = xmlobject + if not (xmlobject is None): + for curve_obj in xmlobject: + for curve_property in curve_obj: + if curve_property.tag == "name": + curve_name = curve_property.text + self.__dict__[curve_name] = curve.Curve(curve_obj) + + def __checkcurve(self, dictionary): + required = ["coords", "values"] + optional = ["name"] + for k, v in dictionary.items(): + if not k in (required+optional): + raise RuntimeError(f"{k} is not a valid property field for a curve.") + for entry in required: + if not entry in dictionary: + raise RuntimeError(f"{entry} is required for creating a curve.") + + def __setitem__(self, key, item): + if not isinstance(item, dict): + raise RuntimeError("Item must be a dictionary") + if len(item) == 0: + self.__delitem__(key) + return + self.__checkcurve(item) + if key in self.__dict__: + self.__delitem__(key) + curve_obj = ET.SubElement(self.xmlobject, "curve") + assert(len(item["coords"])==len(item["values"])) + for i, entry in enumerate(item["coords"]): + item["coords"][i] = str(entry) + item["values"][i] = str(item["values"][i]) + q = ET.SubElement(curve_obj, "name") + q.text = key + q = ET.SubElement(curve_obj, "coords") + q.text = ' '.join(item["coords"]) + q = ET.SubElement(curve_obj, "values") + q.text = ' '.join(item["values"]) + return curve_obj + + + def __getitem__(self, key): + if not (key in ["tree", "xmlobject"]): + return self.__dict__[key] + + def __repr__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree", "name", "xmlobject"]): + newdict[k] = v + return repr(newdict) + + def __len__(self): + return len(self.__dict__) + + def __delitem__(self, key): + obj = self.__dict__[key].xmlobject + obj.getparent().remove(obj) + del self.__dict__[key] + + def clear(self): + return self.__dict__.clear() + + def copy(self): + return self.__dict__.copy() + + def has_key(self, k): + if not (k in ["tree","xmlobject"]): + return k in self.__dict__ + + def update(self, *args, **kwargs): + pass + # return self.__dict__.update(*args, **kwargs) + + def keys(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree","xmlobject"]): + newdict[k] = v + return newdict.keys() + + def values(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree", "xmlobject"]): + newdict[k] = v + return newdict.values() + + def items(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree", "xmlobject"]): + newdict[k] = v + return newdict.items() + + def pop(self, *args): + pass + #return self.__dict__.pop(*args) + + def __cmp__(self, dict_): + return self.__cmp__(self.__dict__, dict_) + + def __contains__(self, item): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree", "xmlobject"]): + newdict[k] = v + return item in newdict + + def __iter__(self): + newdict = {} + for k, v in self.__dict__.items(): + if not (k in ["tree", "xmlobject"]): + newdict[k] = v + return iter(newdict) def add_curve(self, **args): """ diff --git a/ogs6py/classes/parameter_type.py b/ogs6py/classes/parameter_type.py index 42f6084..cb48ed9 100644 --- a/ogs6py/classes/parameter_type.py +++ b/ogs6py/classes/parameter_type.py @@ -7,8 +7,18 @@ """ from fastcore.utils import * +import numpy as np from lxml import etree as ET - +try: + import vtuIO + has_vtuinterface = True +except ImportError: + has_vtuinterface = False +try: + import cexprtk + has_cexprtk = True +except ImportError: + has_cexprtk = False # pylint: disable=C0103, R0902, R0914, R0913 @@ -38,7 +48,7 @@ def __setitem__(self, key, item): expression_counter += 1 parameter_property.text = str(item[expression_counter]) else: - parameter_property.text = (item) + parameter_property.text = str(item) def __getitem__(self, key): if not (key == "xmlobject"): @@ -115,24 +125,91 @@ def __iter__(self): return iter(newdict) class Constant(Parameter_type): - pass + def evaluate_values(self): + if "values" in self.__dict__: + return self.__dict__["values"] + else: + return self.__dict__["value"] class Function(Parameter_type): - pass + def evaluate_values(self): + if has_vtuinterface is False: + raise RuntimeError("VTUinterface is not installed") + if has_cexprtk is False: + raise RuntimeError("cexprtk is not installed") + try: + mesh = self.__dict__["mesh"] + meshfiles = self.xmlobject.getparent().getparent().findall("./mesh") + for file in meshfiles: + if mesh in file.text: + meshfile = file.text + except KeyError: + meshfile = self.xmlobject.getparent().getparent().find("./mesh") + m = vtuIO.VTUIO(meshfile) + st = cexprtk.Symbol_Table({'x': 0.0, 'y': 0.0, 'z': 0.0}, add_constants=True) + dim1 = len(m.points) + dim2 = len(self.__dict__["expression"]) + if dim2 == 1: + array = np.zeros(dim1) + evaluate = cexprtk.Expression(self.__dict__["expression"][0], st) + for i in range(dim1): + st.variables["x"] = m.points[i][0] + st.variables["y"] = m.points[i][1] + st.variables["z"] = m.points[i][2] + array[i] = evaluate() + else: + array = np.zeros((dim1, dim2)) + evaluate = [] + for i in range(dim2): + evaluate.append(cexprtk.Expression(self.__dict__["expression"][i], st)) + for i in range(dim1): + for j in range(dim2): + st.variables["x"] = m.points[i][0] + st.variables["y"] = m.points[i][1] + st.variables["z"] = m.points[i][2] + array[i,j]=evaluate[j]() + return m.points, array, m + class MeshNode(Parameter_type): - pass + def evaluate_values(self): + if has_vtuinterface is False: + raise RuntimeError("vtuIO is not installed") + try: + mesh = self.__dict__["mesh"] + meshfiles = self.xmlobject.getparent().getparent().findall("./mesh") + for file in meshfiles: + if mesh in file.text: + meshfile = file.text + except KeyError: + meshfile = self.xmlobject.getparent().getparent().find("./mesh") + m = vtuIO.VTUIO(meshfile) + array = m.get_point_field(self.__dict__["field_name"]) + return m.points, array, m class MeshElement(Parameter_type): - pass + def evaluate_values(self): + if has_vtuinterface is False: + raise RuntimeError("vtuIO is not installed") + try: + mesh = self.__dict__["mesh"] + meshfiles = self.xmlobject.getparent().getparent().findall("./mesh") + for file in meshfiles: + if mesh in file.text: + meshfile = file.text + except KeyError: + meshfile = self.xmlobject.getparent().getparent().find("./mesh") + m = vtuIO.VTUIO(meshfile) + array = m.get_cell_field(self.__dict__["field_name"]) + return m.cell_center_points, array, m class CurveScaled(Parameter_type): pass -class TimeDependentHeterogeneousParameter(Parameter_type): - pass +#class TimeDependentHeterogeneousParameter(Parameter_type): +# pass class RandomFieldMeshElementParameter(Parameter_type): pass -class Group(Parameter_type): - pass \ No newline at end of file +#class Group(Parameter_type): +# pass \ No newline at end of file diff --git a/ogs6py/classes/parameters.py b/ogs6py/classes/parameters.py index d9664eb..d6aa772 100644 --- a/ogs6py/classes/parameters.py +++ b/ogs6py/classes/parameters.py @@ -7,6 +7,7 @@ """ # pylint: disable=C0103, R0902, R0914, R0913 +from lxml import etree as ET from ogs6py.classes import build_tree from ogs6py.classes import parameter_type @@ -25,7 +26,6 @@ def __init__(self, xmlobject=None): } self.parameter = {} self.xmlobject = xmlobject - #self.__dict__ = {} if not (xmlobject is None): for prmt in xmlobject: for parameter_property in prmt: @@ -43,12 +43,78 @@ def __init__(self, xmlobject=None): self.__dict__[param_name] = parameter_type.MeshElement(prmt) elif param_type == "CurveScaled": self.__dict__[param_name] = parameter_type.CurveScaled(prmt) - elif param_type == "TimeDependentHeterogeneousParameter": - self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) +# elif param_type == "TimeDependentHeterogeneousParameter": +# self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) elif param_type == "RandomFieldMeshElementParameter": self.__dict__[param_name] = parameter_type.RandomFieldMeshElementParameter(prmt) - elif param_type == "Group": - self.__dict__[param_name] = parameter_type.Group(prmt) +# elif param_type == "Group": +# self.__dict__[param_name] = parameter_type.Group(prmt) + + def __checkparameter(self, dictionary): + required = {"Constant": ["name", "type"], + "Function": ["name", "type", "expression"], + "MeshNode": ["name", "type", "field_name"], + "MeshElement": ["name", "type", "field_name"], + "CurvedScaled": ["name", "type", "curve", "parameter"], + "TimeDependentHeterogeneousParameter": ["name", "type", "time_series"], + "RandomFieldMeshElementParameter": ["name", "type","field_name", "range", "seed"], + "Group": ["name", "type", "group_id_property"]} + optional = {"Constant": ["value", "values"], + "Function": ["mesh"], + "MeshNode": ["mesh"], + "MeshElement": ["mesh"], + "CurvedScaled": ["mesh"], + "TimeDependentHeterogeneousParameter": ["mesh"], + "RandomFieldMeshElementParameter": ["mesh"], + "Group": ["mesh"]} + for k, v in dictionary.items(): + if not k in (required[dictionary["type"]]+optional[dictionary["type"]]): + raise RuntimeError(f"{k} is not a valid property field for the specified type.") + for entry in required[dictionary["type"]]: + if not entry in dictionary: + raise RuntimeError(f"{entry} is required for the specified type.") + if dictionary["type"] == "Constant": + if not (("value" in dictionary) or ("values" in dictionary)): + raise RuntimeError("The Constant parameter requires value or values to be specified.") + + + def __setitem__(self, key, item): + if not isinstance(item, dict): + raise RuntimeError("Item must be a dictionary") + if len(item) == 0: + self.__delitem__(key) + return + self.__checkparameter(item) + if key in self.__dict__: + self.__delitem__(key) + prmt_obj = ET.SubElement(self.xmlobject, "parameter") + for k, v in item.items(): + if k == "expression": + q = [] + for subentry in v: + q.append(ET.SubElement(prmt_obj, "expression")) + q[-1].text = subentry + else: + q = ET.SubElement(prmt_obj, k) + q.text = v + if item["type"] == "Constant": + self.__dict__[key] = parameter_type.Constant(prmt_obj) + elif item["type"] == "Function": + self.__dict__[key] = parameter_type.Function(prmt_obj) + elif item["type"] == "MeshNode": + self.__dict__[key] = parameter_type.MeshNode(prmt_obj) + elif item["type"] == "MeshElement": + self.__dict__[key] = parameter_type.MeshElement(prmt_obj) + elif item["type"] == "CurveScaled": + self.__dict__[key] = parameter_type.CurveScaled(prmt_obj) +# elif item["type"] == "TimeDependentHeterogeneousParameter": +# self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) + elif item["type"] == "RandomFieldMeshElementParameter": + self.__dict__[key] = parameter_type.RandomFieldMeshElementParameter(prmt_obj) +# elif item["type"] == "Group": +# self.__dict__[param_name] = parameter_type.Group(prmt) + return prmt_obj + def __getitem__(self, key): if not (key in ["tree","parameter", "xmlobject"]): @@ -57,7 +123,7 @@ def __getitem__(self, key): def __repr__(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "xmlobject"]): + if not (k in ["tree","parameter", "name", "xmlobject"]): newdict[k] = v return repr(newdict) @@ -65,8 +131,9 @@ def __len__(self): return len(self.__dict__) def __delitem__(self, key): - pass - #del self.__dict__[key] + obj = self.__dict__[key].xmlobject + obj.getparent().remove(obj) + del self.__dict__[key] def clear(self): return self.__dict__.clear() diff --git a/ogs6py/ogs.py b/ogs6py/ogs.py index 90d5306..578cd91 100644 --- a/ogs6py/ogs.py +++ b/ogs6py/ogs.py @@ -50,8 +50,8 @@ def __init__(self, input_file=None, project_file="default.prj", xmlstring=None, self.media = media.Media() self.timeloop = timeloop.TimeLoop() self.local_coordinate_system = local_coordinate_system.LocalCoordinateSystem() - #self.parameters = parameters.Parameters() - self.curves = curves.Curves() + self.__parameters = None + self.__curves = None self.processvars = processvars.ProcessVars() self.linsolvers = linsolvers.LinSolvers() self.nonlinsolvers = nonlinsolvers.NonLinSolvers() @@ -74,7 +74,7 @@ def __init__(self, input_file=None, project_file="default.prj", xmlstring=None, self.omp_num_threads = None if "PROJECT_FILE" in args: self.project_file = args['PROJECT_FILE'] - else: + if self.project_file is None: print("PROJECT_FILE for output not given. Calling it default.prj.") self.project_file = "default.prj" if "INPUT_FILE" in args: @@ -86,14 +86,8 @@ def __init__(self, input_file=None, project_file="default.prj", xmlstring=None, display.Display(self.tree) else: raise RuntimeError(f"Input project file {args['INPUT_FILE']} not found.") - if "XMLSTRING" in args or (self.xmlstring is not None): - root = ET.fromstring(args['XMLSTRING']) - self.tree = ET.ElementTree(root) - try: - paramobj = self.tree.find("./parameters") - except AttributeError: - paramobj = None - self.parameters = parameters.Parameters(xmlobject=paramobj) + if "XMLSTRING" in args: + self.xmlstring = args["XMLSTRING"] def __dict2xml(self, parent, dictionary): for entry in dictionary: @@ -124,6 +118,9 @@ def _get_root(self): if self.tree is None: if self.input_file is not None: self.tree = ET.parse(self.input_file) + elif self.xmlstring is not None: + root = ET.fromstring(self.xmlstring) + self.tree = ET.ElementTree(root) else: self.build_tree() root = self.tree.getroot() @@ -143,6 +140,26 @@ def _get_root(self): self.include_elements.append(child) return root + @property + def parameters(self): + try: + paramobj = self.tree.find("./parameters") + except AttributeError: + paramobj = None + if not (paramobj == self.__parameters): + self.__parameters = parameters.Parameters(xmlobject=paramobj) + return self.__parameters + + @property + def curves(self): + try: + curveobj = self.tree.find("./curves") + except AttributeError: + curveobj = None + if not (curveobj == self.__curves): + self.__curves = curves.Curves(xmlobject=curveobj) + return self.__curves + @classmethod def _get_parameter_pointer(cls, root, name, xpath): params = root.findall(xpath) diff --git a/setup.py b/setup.py index 1ea70ae..0349813 100644 --- a/setup.py +++ b/setup.py @@ -55,5 +55,6 @@ def find_version(*file_paths): include_package_data=True, python_requires='>=3.8', install_requires=["lxml","pandas"], + extras_require={"parameter_access": ["VTUinterface", "cexprtk"]}, py_modules=["ogs6py/ogs","ogs6py/log_parser/log_parser", "ogs6py/log_parser/common_ogs_analyses", "ogs6py/ogs_regexes/ogs_regexes"], packages=["ogs6py/classes","ogs6py/log_parser","ogs6py/ogs_regexes"]) From dcc2e3e24a42f70a48856f90547d51edc1eb27af Mon Sep 17 00:00:00 2001 From: joergbuchwald Date: Tue, 7 Feb 2023 22:22:59 +0100 Subject: [PATCH 3/6] fix from-scratch-add_parameter issue --- ogs6py/classes/parameters.py | 35 +++++++++++++++++++++++++++-------- ogs6py/ogs.py | 12 ++++++++---- 2 files changed, 35 insertions(+), 12 deletions(-) diff --git a/ogs6py/classes/parameters.py b/ogs6py/classes/parameters.py index d6aa772..d7dcb32 100644 --- a/ogs6py/classes/parameters.py +++ b/ogs6py/classes/parameters.py @@ -16,6 +16,7 @@ class Parameters(build_tree.BuildTree): Class for managing the parameters section of the project file. """ def __init__(self, xmlobject=None): + #print("init") self.tree = { 'parameters': { 'tag': 'parameters', @@ -51,6 +52,7 @@ def __init__(self, xmlobject=None): # self.__dict__[param_name] = parameter_type.Group(prmt) def __checkparameter(self, dictionary): + #print("checkparameter") required = {"Constant": ["name", "type"], "Function": ["name", "type", "expression"], "MeshNode": ["name", "type", "field_name"], @@ -79,6 +81,7 @@ def __checkparameter(self, dictionary): def __setitem__(self, key, item): + #print("setitem") if not isinstance(item, dict): raise RuntimeError("Item must be a dictionary") if len(item) == 0: @@ -117,77 +120,92 @@ def __setitem__(self, key, item): def __getitem__(self, key): - if not (key in ["tree","parameter", "xmlobject"]): + #print("getitem") + if not (key in ["tree", "parameter", "xmlobject"]): return self.__dict__[key] def __repr__(self): + #print("repr") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "name", "xmlobject"]): + if not (k in ["tree", "parameter", "name", "xmlobject"]): newdict[k] = v return repr(newdict) def __len__(self): + #print("len") return len(self.__dict__) def __delitem__(self, key): + #print("delitem") obj = self.__dict__[key].xmlobject obj.getparent().remove(obj) del self.__dict__[key] def clear(self): + #print("clear") return self.__dict__.clear() def copy(self): + #print("copy") return self.__dict__.copy() def has_key(self, k): - if not (k in ["tree","parameter", "xmlobject"]): + #print("has key") + if not (k in ["tree", "parameter", "xmlobject"]): return k in self.__dict__ def update(self, *args, **kwargs): + #print("update") pass # return self.__dict__.update(*args, **kwargs) def keys(self): + #print("keys") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject"]): newdict[k] = v return newdict.keys() def values(self): + #print("values") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject"]): newdict[k] = v return newdict.values() def items(self): + #print("items") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject"]): newdict[k] = v return newdict.items() def pop(self, *args): + #print("pop") pass #return self.__dict__.pop(*args) def __cmp__(self, dict_): + #print("cmp") return self.__cmp__(self.__dict__, dict_) def __contains__(self, item): + #print("contains") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject"]): newdict[k] = v return item in newdict def __iter__(self): + #print("iter") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree","parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject"]): newdict[k] = v return iter(newdict) @@ -210,6 +228,7 @@ def add_parameter(self, **args): parameter_name : `list` use_local_coordinate_system : `bool` or `str` """ + #print("add parameter") self._convertargs(args) if "name" not in args: raise KeyError("No parameter name given.") diff --git a/ogs6py/ogs.py b/ogs6py/ogs.py index 578cd91..1d3c169 100644 --- a/ogs6py/ogs.py +++ b/ogs6py/ogs.py @@ -144,20 +144,24 @@ def _get_root(self): def parameters(self): try: paramobj = self.tree.find("./parameters") + if not (paramobj == self.__parameters): + self.__parameters = parameters.Parameters(xmlobject=paramobj) except AttributeError: paramobj = None - if not (paramobj == self.__parameters): - self.__parameters = parameters.Parameters(xmlobject=paramobj) + if self.__parameters is None: + self.__parameters = parameters.Parameters(xmlobject=paramobj) return self.__parameters @property def curves(self): try: curveobj = self.tree.find("./curves") + if not (curveobj == self.__curves): + self.__curves = curves.Curves(xmlobject=curveobj) except AttributeError: curveobj = None - if not (curveobj == self.__curves): - self.__curves = curves.Curves(xmlobject=curveobj) + if self.__curves is None: + self.__curves = curves.Curves(xmlobject=curveobj) return self.__curves @classmethod From a0246944a63819496cbe92d70b42c359e5340439 Mon Sep 17 00:00:00 2001 From: joergbuchwald Date: Thu, 9 Feb 2023 18:12:35 +0100 Subject: [PATCH 4/6] implement curves in Function parameter --- ogs6py/classes/parameter_type.py | 37 ++++++++++++------- ogs6py/classes/parameters.py | 62 ++++++++++++-------------------- ogs6py/ogs.py | 4 +-- 3 files changed, 48 insertions(+), 55 deletions(-) diff --git a/ogs6py/classes/parameter_type.py b/ogs6py/classes/parameter_type.py index cb48ed9..b0e7e94 100644 --- a/ogs6py/classes/parameter_type.py +++ b/ogs6py/classes/parameter_type.py @@ -23,9 +23,10 @@ # pylint: disable=C0103, R0902, R0914, R0913 class Parameter_type: - def __init__(self, xmlobject:object=None) -> None: + def __init__(self, xmlobject:object=None, curvesobj=None) -> None: self.__dict__ = {} self.xmlobject = xmlobject + self.curvesobj = curvesobj if not self.xmlobject is None: for parameter_property in self.xmlobject: if parameter_property.tag == "expression": @@ -51,13 +52,13 @@ def __setitem__(self, key, item): parameter_property.text = str(item) def __getitem__(self, key): - if not (key == "xmlobject"): + if not (key in ["xmlobject", "curvesobj"]): return self.__dict__[key] def __repr__(self): newdict = {} for k, v in self.__dict__.items(): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): newdict[k] = v return repr(newdict) @@ -75,7 +76,7 @@ def copy(self): return self.__dict__.copy() def has_key(self, k): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): return k in self.__dict__ def update(self, *args, **kwargs): @@ -85,21 +86,21 @@ def update(self, *args, **kwargs): def keys(self): newdict = {} for k, v in self.__dict__.items(): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): newdict[k] = v return newdict.keys() def values(self): newdict = {} for k, v in self.__dict__.items(): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): newdict[k] = v return newdict.values() def items(self): newdict = {} for k, v in self.__dict__.items(): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): newdict[k] = v return newdict.items() @@ -113,14 +114,14 @@ def __cmp__(self, dict_): def __contains__(self, item): newdict = {} for k, v in self.__dict__.items(): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): newdict[k] = v return item in newdict def __iter__(self): newdict = {} for k, v in self.__dict__.items(): - if not (k == "xmlobject"): + if not (k in ["xmlobject", "curvesobj"]): newdict[k] = v return iter(newdict) @@ -132,7 +133,7 @@ def evaluate_values(self): return self.__dict__["value"] class Function(Parameter_type): - def evaluate_values(self): + def evaluate_values(self, t=0): if has_vtuinterface is False: raise RuntimeError("VTUinterface is not installed") if has_cexprtk is False: @@ -144,9 +145,14 @@ def evaluate_values(self): if mesh in file.text: meshfile = file.text except KeyError: - meshfile = self.xmlobject.getparent().getparent().find("./mesh") + meshfile = self.xmlobject.getparent().getparent().find("./mesh").text m = vtuIO.VTUIO(meshfile) - st = cexprtk.Symbol_Table({'x': 0.0, 'y': 0.0, 'z': 0.0}, add_constants=True) + st = cexprtk.Symbol_Table({'x': 0.0, 'y': 0.0, 'z': 0.0, 't': t}, add_constants=True) + try: + for curve in self.curvesobj.keys(): + st.functions[curve] = self.curvesobj[curve].evaluate_values + except: + pass dim1 = len(m.points) dim2 = len(self.__dict__["expression"]) if dim2 == 1: @@ -204,7 +210,12 @@ def evaluate_values(self): return m.cell_center_points, array, m class CurveScaled(Parameter_type): - pass + def evaluate_values(self, curve_coords=None): + if curve_coords is None: + t_start = self.xmlobject.getparent().getparent().find("./time_loop/processes/process/time_stepping/t_initial") + t_end = self.xmlobject.getparent().getparent().find("./time_loop/processes/process/time_stepping/t_end") + curve_coords = np.linspace(t_start, t_end, 1000, endpoint=True) + return #class TimeDependentHeterogeneousParameter(Parameter_type): # pass diff --git a/ogs6py/classes/parameters.py b/ogs6py/classes/parameters.py index d7dcb32..9c88772 100644 --- a/ogs6py/classes/parameters.py +++ b/ogs6py/classes/parameters.py @@ -15,8 +15,7 @@ class Parameters(build_tree.BuildTree): """ Class for managing the parameters section of the project file. """ - def __init__(self, xmlobject=None): - #print("init") + def __init__(self, xmlobject=None, curvesobj=None): self.tree = { 'parameters': { 'tag': 'parameters', @@ -27,6 +26,7 @@ def __init__(self, xmlobject=None): } self.parameter = {} self.xmlobject = xmlobject + self.curvesobj = curvesobj if not (xmlobject is None): for prmt in xmlobject: for parameter_property in prmt: @@ -35,24 +35,23 @@ def __init__(self, xmlobject=None): elif parameter_property.tag == "name": param_name = parameter_property.text if param_type == "Constant": - self.__dict__[param_name] = parameter_type.Constant(prmt) + self.__dict__[param_name] = parameter_type.Constant(prmt, curvesobj) elif param_type == "Function": - self.__dict__[param_name] = parameter_type.Function(prmt) + self.__dict__[param_name] = parameter_type.Function(prmt, curvesobj) elif param_type == "MeshNode": - self.__dict__[param_name] = parameter_type.MeshNode(prmt) + self.__dict__[param_name] = parameter_type.MeshNode(prmt, curvesobj) elif param_type == "MeshElement": - self.__dict__[param_name] = parameter_type.MeshElement(prmt) + self.__dict__[param_name] = parameter_type.MeshElement(prmt, curvesobj) elif param_type == "CurveScaled": - self.__dict__[param_name] = parameter_type.CurveScaled(prmt) + self.__dict__[param_name] = parameter_type.CurveScaled(prmt, curvesobj) # elif param_type == "TimeDependentHeterogeneousParameter": # self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) elif param_type == "RandomFieldMeshElementParameter": - self.__dict__[param_name] = parameter_type.RandomFieldMeshElementParameter(prmt) + self.__dict__[param_name] = parameter_type.RandomFieldMeshElementParameter(prmt, curvesobj) # elif param_type == "Group": # self.__dict__[param_name] = parameter_type.Group(prmt) def __checkparameter(self, dictionary): - #print("checkparameter") required = {"Constant": ["name", "type"], "Function": ["name", "type", "expression"], "MeshNode": ["name", "type", "field_name"], @@ -81,7 +80,6 @@ def __checkparameter(self, dictionary): def __setitem__(self, key, item): - #print("setitem") if not isinstance(item, dict): raise RuntimeError("Item must be a dictionary") if len(item) == 0: @@ -101,111 +99,96 @@ def __setitem__(self, key, item): q = ET.SubElement(prmt_obj, k) q.text = v if item["type"] == "Constant": - self.__dict__[key] = parameter_type.Constant(prmt_obj) + self.__dict__[key] = parameter_type.Constant(prmt_obj, self.curvesobj) elif item["type"] == "Function": - self.__dict__[key] = parameter_type.Function(prmt_obj) + self.__dict__[key] = parameter_type.Function(prmt_obj, self.curvesobj) elif item["type"] == "MeshNode": - self.__dict__[key] = parameter_type.MeshNode(prmt_obj) + self.__dict__[key] = parameter_type.MeshNode(prmt_obj, self.curvesobj) elif item["type"] == "MeshElement": - self.__dict__[key] = parameter_type.MeshElement(prmt_obj) + self.__dict__[key] = parameter_type.MeshElement(prmt_obj, self.curvesobj) elif item["type"] == "CurveScaled": - self.__dict__[key] = parameter_type.CurveScaled(prmt_obj) + self.__dict__[key] = parameter_type.CurveScaled(prmt_obj, self.curvesobj) # elif item["type"] == "TimeDependentHeterogeneousParameter": # self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) elif item["type"] == "RandomFieldMeshElementParameter": - self.__dict__[key] = parameter_type.RandomFieldMeshElementParameter(prmt_obj) + self.__dict__[key] = parameter_type.RandomFieldMeshElementParameter(prmt_obj, self.curvesobj) # elif item["type"] == "Group": # self.__dict__[param_name] = parameter_type.Group(prmt) return prmt_obj def __getitem__(self, key): - #print("getitem") - if not (key in ["tree", "parameter", "xmlobject"]): + if not (key in ["tree", "parameter", "xmlobject", "curvesobj"]): return self.__dict__[key] def __repr__(self): - #print("repr") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "name", "xmlobject"]): + if not (k in ["tree", "parameter", "name", "xmlobject", "curvesobj"]): newdict[k] = v return repr(newdict) def __len__(self): - #print("len") return len(self.__dict__) def __delitem__(self, key): - #print("delitem") obj = self.__dict__[key].xmlobject obj.getparent().remove(obj) del self.__dict__[key] def clear(self): - #print("clear") return self.__dict__.clear() def copy(self): - #print("copy") return self.__dict__.copy() def has_key(self, k): - #print("has key") - if not (k in ["tree", "parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): return k in self.__dict__ def update(self, *args, **kwargs): - #print("update") pass # return self.__dict__.update(*args, **kwargs) def keys(self): - #print("keys") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): newdict[k] = v return newdict.keys() def values(self): - #print("values") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): newdict[k] = v return newdict.values() def items(self): - #print("items") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): newdict[k] = v return newdict.items() def pop(self, *args): - #print("pop") pass #return self.__dict__.pop(*args) def __cmp__(self, dict_): - #print("cmp") return self.__cmp__(self.__dict__, dict_) def __contains__(self, item): - #print("contains") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): newdict[k] = v return item in newdict def __iter__(self): - #print("iter") newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): newdict[k] = v return iter(newdict) @@ -228,7 +211,6 @@ def add_parameter(self, **args): parameter_name : `list` use_local_coordinate_system : `bool` or `str` """ - #print("add parameter") self._convertargs(args) if "name" not in args: raise KeyError("No parameter name given.") diff --git a/ogs6py/ogs.py b/ogs6py/ogs.py index 1d3c169..789f581 100644 --- a/ogs6py/ogs.py +++ b/ogs6py/ogs.py @@ -145,11 +145,11 @@ def parameters(self): try: paramobj = self.tree.find("./parameters") if not (paramobj == self.__parameters): - self.__parameters = parameters.Parameters(xmlobject=paramobj) + self.__parameters = parameters.Parameters(xmlobject=paramobj, curvesobj=self.curves) except AttributeError: paramobj = None if self.__parameters is None: - self.__parameters = parameters.Parameters(xmlobject=paramobj) + self.__parameters = parameters.Parameters(xmlobject=paramobj, curvesobj=self.curves) return self.__parameters @property From bfe46c976fff52063d57a84fdb9f6cdd36d31b1c Mon Sep 17 00:00:00 2001 From: joergbuchwald Date: Fri, 10 Feb 2023 16:13:35 +0100 Subject: [PATCH 5/6] implementation of CurveScaled parameter --- ogs6py/classes/parameter_type.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/ogs6py/classes/parameter_type.py b/ogs6py/classes/parameter_type.py index b0e7e94..7eedf26 100644 --- a/ogs6py/classes/parameter_type.py +++ b/ogs6py/classes/parameter_type.py @@ -35,6 +35,11 @@ def __init__(self, xmlobject:object=None, curvesobj=None) -> None: except: self.__dict__["expression"] = [] self.__dict__["expression"].append(parameter_property.text) + elif parameter_property.tag == "value": + self.__dict__[parameter_property.tag] = float(parameter_property.text) + elif parameter_property.tag == "values": + # convert to tensor + include local coordinate system + self.__dict__[parameter_property.tag] = np.fromstring(parameter_property.text, sep=' ') else: self.__dict__[parameter_property.tag] = parameter_property.text def __setitem__(self, key, item): @@ -215,7 +220,17 @@ def evaluate_values(self, curve_coords=None): t_start = self.xmlobject.getparent().getparent().find("./time_loop/processes/process/time_stepping/t_initial") t_end = self.xmlobject.getparent().getparent().find("./time_loop/processes/process/time_stepping/t_end") curve_coords = np.linspace(t_start, t_end, 1000, endpoint=True) - return + parameter_name = self.__dict__["parameter"] + curve_name = self.__dict__["curve"] + parameter_type = self.xmlobject.getparent().find(f"./parameter[name='{parameter_name}'/type").text + try: + parameter_value = float(self.xmlobject.getparent().find(f"./parameter[name='{parameter_name}'/value").text) + except: + raise RuntimeError("This function is implemented only for constant single parameter values") + if parameter_type == "Constant": + return parameter_value*self.curvesobj[curve_name].evaluate_values(curve_coords) + else: + raise RuntimeError("This function is implemented constant parameter types only") #class TimeDependentHeterogeneousParameter(Parameter_type): # pass From 63af8adca8fb0ebaf7f2e3f2c398d30dce1526bf Mon Sep 17 00:00:00 2001 From: joergbuchwald Date: Tue, 14 Feb 2023 22:22:38 +0100 Subject: [PATCH 6/6] implement local coordinate system --- ogs6py/classes/local_coordinate_system.py | 18 +++++++- ogs6py/classes/parameter_type.py | 54 +++++++++++++++++------ ogs6py/classes/parameters.py | 46 ++++++++++--------- ogs6py/ogs.py | 48 +++++++++++++------- 4 files changed, 115 insertions(+), 51 deletions(-) diff --git a/ogs6py/classes/local_coordinate_system.py b/ogs6py/classes/local_coordinate_system.py index e82e15d..da56a7b 100644 --- a/ogs6py/classes/local_coordinate_system.py +++ b/ogs6py/classes/local_coordinate_system.py @@ -6,6 +6,7 @@ http://www.opengeosys.org/project/license """ +import numpy as np # pylint: disable=C0103, R0902, R0914, R0913 from ogs6py.classes import build_tree @@ -13,7 +14,7 @@ class LocalCoordinateSystem(build_tree.BuildTree): """ Class for defining a local coordinate system in the project file" """ - def __init__(self): + def __init__(self, xmlobject=None): self.tree = { 'local_coordinate_system': { 'tag': 'local_coordinate_system', @@ -22,6 +23,21 @@ def __init__(self): 'children': {} } } + self.xmlobject = xmlobject + self.R = None + if not self.xmlobject is None: + basis_vectors = self.xmlobject.getchildren() + dim = len(basis_vectors) + self.R = np.zeros((dim,dim)) + basis_vector_names = [] + for vec in basis_vectors: + basis_vector_names.append(vec.text) + basis_vector_values = [] + for vec in basis_vector_names: + basis_vector_values.append(np.fromstring(self.xmlobject.getparent().find(f"./parameters/parameter[name='{vec}']/values").text, sep=' ')) + for i in range(dim): + for j in range(dim): + self.R[i,j] = basis_vector_values[i][j] def add_basis_vec(self, **args): """ diff --git a/ogs6py/classes/parameter_type.py b/ogs6py/classes/parameter_type.py index 7eedf26..9a6d469 100644 --- a/ogs6py/classes/parameter_type.py +++ b/ogs6py/classes/parameter_type.py @@ -8,6 +8,7 @@ """ from fastcore.utils import * import numpy as np +import ctypes from lxml import etree as ET try: import vtuIO @@ -23,10 +24,12 @@ # pylint: disable=C0103, R0902, R0914, R0913 class Parameter_type: - def __init__(self, xmlobject:object=None, curvesobj=None) -> None: + def __init__(self, xmlobject:object=None, paramobject=None, curvesobj=None, trafo_matrix=None) -> None: self.__dict__ = {} self.xmlobject = xmlobject + self.paramobject = id(paramobject) self.curvesobj = curvesobj + self.trafo_matrix = trafo_matrix if not self.xmlobject is None: for parameter_property in self.xmlobject: if parameter_property.tag == "expression": @@ -57,13 +60,13 @@ def __setitem__(self, key, item): parameter_property.text = str(item) def __getitem__(self, key): - if not (key in ["xmlobject", "curvesobj"]): + if not (key in ["xmlobject", "curvesobj", "trafo_matrix"]): return self.__dict__[key] def __repr__(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return repr(newdict) @@ -81,7 +84,7 @@ def copy(self): return self.__dict__.copy() def has_key(self, k): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix", "paramobject"]): return k in self.__dict__ def update(self, *args, **kwargs): @@ -91,21 +94,21 @@ def update(self, *args, **kwargs): def keys(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix", "paramobject"]): newdict[k] = v return newdict.keys() def values(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix", "paramobject"]): newdict[k] = v return newdict.values() def items(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix", "paramobject"]): newdict[k] = v return newdict.items() @@ -119,21 +122,42 @@ def __cmp__(self, dict_): def __contains__(self, item): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix", "paramobject"]): newdict[k] = v return item in newdict def __iter__(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["xmlobject", "curvesobj"]): + if not (k in ["xmlobject", "curvesobj", "trafo_matrix", "paramobject"]): newdict[k] = v return iter(newdict) class Constant(Parameter_type): def evaluate_values(self): if "values" in self.__dict__: - return self.__dict__["values"] + values_size = len(self.__dict__["values"]) + if values_size == 1: + return self.__dict__["values"] + elif "local_coordinate_system" in self.__dict__: + if self.__dict__["local_coordinate_system"] == "true": + dim_trafo = self.trafo_matrix.shape[0] + if values_size == dim_trafo: + return np.matmul(self.trafo_matrix, self.__dict__["values"]) + elif values_size == dim_trafo**2: + values_matrix = np.zeros((dim_trafo,dim_trafo)) + mapping = {} + mapping[2] = {0: (0,0), 1: (0,1), 2: (1,0), 3: (1,1)} + mapping[3] = {0: (0,0), 1: (0,1), 2: (0,2), + 3: (1,0), 4: (1,1), 5: (1,2), + 6: (2,0), 7: (2,1), 8: (2,2)} + for i, val in enumerate(self.__dict__["values"]): + values_matrix[mapping[dim_trafo][i][0],mapping[dim_trafo][i][1]] = val + return np.matmul(self.trafo_matrix,np.matmul(values_matrix,self.trafo_matrix.transpose())) + else: + raise RuntimeError("Parameter size in combination with local coordinate system is not supported yet.") + else: + return self.__dict__["values"] else: return self.__dict__["value"] @@ -224,11 +248,15 @@ def evaluate_values(self, curve_coords=None): curve_name = self.__dict__["curve"] parameter_type = self.xmlobject.getparent().find(f"./parameter[name='{parameter_name}'/type").text try: - parameter_value = float(self.xmlobject.getparent().find(f"./parameter[name='{parameter_name}'/value").text) + paramobject = ctypes.cast(self.paramobject, ctypes.py_object).value + parameter_value = paramobject[parameter_name].evaluate_values() except: - raise RuntimeError("This function is implemented only for constant single parameter values") + raise RuntimeError("Can't find parameter.") if parameter_type == "Constant": - return parameter_value*self.curvesobj[curve_name].evaluate_values(curve_coords) + if len(parameter_value) == 1: + return parameter_value*self.curvesobj[curve_name].evaluate_values(curve_coords) + else: + return [val*self.curvesobj[curve_name].evaluate_values(curve_coords) for val in parameter_value] else: raise RuntimeError("This function is implemented constant parameter types only") diff --git a/ogs6py/classes/parameters.py b/ogs6py/classes/parameters.py index 9c88772..b87bf4d 100644 --- a/ogs6py/classes/parameters.py +++ b/ogs6py/classes/parameters.py @@ -7,6 +7,7 @@ """ # pylint: disable=C0103, R0902, R0914, R0913 +import numpy as np from lxml import etree as ET from ogs6py.classes import build_tree from ogs6py.classes import parameter_type @@ -15,7 +16,7 @@ class Parameters(build_tree.BuildTree): """ Class for managing the parameters section of the project file. """ - def __init__(self, xmlobject=None, curvesobj=None): + def __init__(self, xmlobject=None, curvesobj=None, trafo_matrix=None): self.tree = { 'parameters': { 'tag': 'parameters', @@ -27,6 +28,7 @@ def __init__(self, xmlobject=None, curvesobj=None): self.parameter = {} self.xmlobject = xmlobject self.curvesobj = curvesobj + self.trafo_matrix = trafo_matrix if not (xmlobject is None): for prmt in xmlobject: for parameter_property in prmt: @@ -35,19 +37,19 @@ def __init__(self, xmlobject=None, curvesobj=None): elif parameter_property.tag == "name": param_name = parameter_property.text if param_type == "Constant": - self.__dict__[param_name] = parameter_type.Constant(prmt, curvesobj) + self.__dict__[param_name] = parameter_type.Constant(prmt, self, curvesobj, trafo_matrix) elif param_type == "Function": - self.__dict__[param_name] = parameter_type.Function(prmt, curvesobj) + self.__dict__[param_name] = parameter_type.Function(prmt, self, curvesobj, trafo_matrix) elif param_type == "MeshNode": - self.__dict__[param_name] = parameter_type.MeshNode(prmt, curvesobj) + self.__dict__[param_name] = parameter_type.MeshNode(prmt, self, curvesobj, trafo_matrix) elif param_type == "MeshElement": - self.__dict__[param_name] = parameter_type.MeshElement(prmt, curvesobj) + self.__dict__[param_name] = parameter_type.MeshElement(prmt, self, curvesobj, trafo_matrix) elif param_type == "CurveScaled": - self.__dict__[param_name] = parameter_type.CurveScaled(prmt, curvesobj) + self.__dict__[param_name] = parameter_type.CurveScaled(prmt, self, curvesobj, trafo_matrix) # elif param_type == "TimeDependentHeterogeneousParameter": # self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) elif param_type == "RandomFieldMeshElementParameter": - self.__dict__[param_name] = parameter_type.RandomFieldMeshElementParameter(prmt, curvesobj) + self.__dict__[param_name] = parameter_type.RandomFieldMeshElementParameter(prmt, self, curvesobj, trafo_matrix) # elif param_type == "Group": # self.__dict__[param_name] = parameter_type.Group(prmt) @@ -99,32 +101,32 @@ def __setitem__(self, key, item): q = ET.SubElement(prmt_obj, k) q.text = v if item["type"] == "Constant": - self.__dict__[key] = parameter_type.Constant(prmt_obj, self.curvesobj) + self.__dict__[key] = parameter_type.Constant(prmt_obj, self, self.curvesobj, self.trafo_matrix) elif item["type"] == "Function": - self.__dict__[key] = parameter_type.Function(prmt_obj, self.curvesobj) + self.__dict__[key] = parameter_type.Function(prmt_obj, self, self.curvesobj, self.trafo_matrix) elif item["type"] == "MeshNode": - self.__dict__[key] = parameter_type.MeshNode(prmt_obj, self.curvesobj) + self.__dict__[key] = parameter_type.MeshNode(prmt_obj, self, self.curvesobj, self.trafo_matrix) elif item["type"] == "MeshElement": - self.__dict__[key] = parameter_type.MeshElement(prmt_obj, self.curvesobj) + self.__dict__[key] = parameter_type.MeshElement(prmt_obj,self, self.curvesobj, self.trafo_matrix) elif item["type"] == "CurveScaled": - self.__dict__[key] = parameter_type.CurveScaled(prmt_obj, self.curvesobj) + self.__dict__[key] = parameter_type.CurveScaled(prmt_obj, self, self.curvesobj, self.trafo_matrix) # elif item["type"] == "TimeDependentHeterogeneousParameter": # self.__dict__[param_name] = parameter_type.TimeDependentHeterogeneousParameter(prmt) elif item["type"] == "RandomFieldMeshElementParameter": - self.__dict__[key] = parameter_type.RandomFieldMeshElementParameter(prmt_obj, self.curvesobj) + self.__dict__[key] = parameter_type.RandomFieldMeshElementParameter(prmt_obj, self, self.curvesobj, self.trafo_matrix) # elif item["type"] == "Group": # self.__dict__[param_name] = parameter_type.Group(prmt) return prmt_obj def __getitem__(self, key): - if not (key in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (key in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): return self.__dict__[key] def __repr__(self): - newdict = {} + newdict = dict() for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "name", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "name", "xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return repr(newdict) @@ -143,7 +145,7 @@ def copy(self): return self.__dict__.copy() def has_key(self, k): - if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): return k in self.__dict__ def update(self, *args, **kwargs): @@ -153,21 +155,21 @@ def update(self, *args, **kwargs): def keys(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return newdict.keys() def values(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return newdict.values() def items(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return newdict.items() @@ -181,14 +183,14 @@ def __cmp__(self, dict_): def __contains__(self, item): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return item in newdict def __iter__(self): newdict = {} for k, v in self.__dict__.items(): - if not (k in ["tree", "parameter", "xmlobject", "curvesobj"]): + if not (k in ["tree", "parameter", "xmlobject", "curvesobj", "trafo_matrix"]): newdict[k] = v return iter(newdict) diff --git a/ogs6py/ogs.py b/ogs6py/ogs.py index 789f581..a6517c6 100644 --- a/ogs6py/ogs.py +++ b/ogs6py/ogs.py @@ -31,7 +31,7 @@ class OGS: Parameters ---------- - project_file : `str`, optional + prjfile : `str`, optional Filename of the output project file Default: default.prj input_file : `str`, optional @@ -42,16 +42,19 @@ class OGS: verbose : `bool`, optional Default: False """ - def __init__(self, input_file=None, project_file="default.prj", xmlstring=None, verbose=False, omp_num_threads=None, **args): + def __init__(self, input_file=None, prjfile="default.prj", xmlstring=None, verbose=False, omp_num_threads=None, **args): self.geo = geo.Geo() self.mesh = mesh.Mesh() self.pyscript = python_script.PythonScript() self.processes = processes.Processes() self.media = media.Media() self.timeloop = timeloop.TimeLoop() - self.local_coordinate_system = local_coordinate_system.LocalCoordinateSystem() + self.__local_coordinate_system = None self.__parameters = None self.__curves = None + self.__local_coordinate_system_obj = None + self.__parameters_obj = None + self.__curves_obj = None self.processvars = processvars.ProcessVars() self.linsolvers = linsolvers.LinSolvers() self.nonlinsolvers = nonlinsolvers.NonLinSolvers() @@ -73,10 +76,10 @@ def __init__(self, input_file=None, project_file="default.prj", xmlstring=None, else: self.omp_num_threads = None if "PROJECT_FILE" in args: - self.project_file = args['PROJECT_FILE'] - if self.project_file is None: + self.prjfile = args['PROJECT_FILE'] + if self.prjfile is None: print("PROJECT_FILE for output not given. Calling it default.prj.") - self.project_file = "default.prj" + self.prjfile = "default.prj" if "INPUT_FILE" in args: self.input_file = args['INPUT_FILE'] if self.input_file is not None: @@ -144,19 +147,21 @@ def _get_root(self): def parameters(self): try: paramobj = self.tree.find("./parameters") - if not (paramobj == self.__parameters): - self.__parameters = parameters.Parameters(xmlobject=paramobj, curvesobj=self.curves) + if not (paramobj == self.__parameters_obj): + self.__parameters_obj = paramobj + self.__parameters = parameters.Parameters(xmlobject=paramobj, curvesobj=self.curves, trafo_matrix=self.local_coordinate_system.R) except AttributeError: paramobj = None if self.__parameters is None: - self.__parameters = parameters.Parameters(xmlobject=paramobj, curvesobj=self.curves) + self.__parameters = parameters.Parameters(xmlobject=paramobj, curvesobj=self.curves, trafo_matrix=self.local_coordinate_system.R) return self.__parameters @property def curves(self): try: curveobj = self.tree.find("./curves") - if not (curveobj == self.__curves): + if not (curveobj == self.__curves_obj): + self.__curves_obj = curveobj self.__curves = curves.Curves(xmlobject=curveobj) except AttributeError: curveobj = None @@ -164,6 +169,19 @@ def curves(self): self.__curves = curves.Curves(xmlobject=curveobj) return self.__curves + @property + def local_coordinate_system(self): + try: + lcsobj = self.tree.find("./local_coordinate_system") + if not (lcsobj == self.__local_coordinate_system_obj): + self.__local_coordinate_system_obj = lcsobj + self.__local_coordinate_system = local_coordinate_system.LocalCoordinateSystem(xmlobject=lcsobj) + except AttributeError: + lcsobj = None + if self.__local_coordinate_system is None: + self.__local_coordinate_system = local_coordinate_system.LocalCoordinateSystem(xmlobject=lcsobj) + return self.__local_coordinate_system + @classmethod def _get_parameter_pointer(cls, root, name, xpath): params = root.findall(xpath) @@ -638,9 +656,9 @@ def run_model(self, logfile="out.log", path=None, args=None, container_path=None if not args is None: cmd += f"{args} " if write_logs is True: - cmd += f"{self.project_file} > {self.logfile}" + cmd += f"{self.prjfile} > {self.logfile}" else: - cmd += f"{self.project_file}" + cmd += f"{self.prjfile}" startt = time.time() if sys.platform == "win32": returncode = subprocess.run(cmd, shell=True, stdout=subprocess.DEVNULL, stderr=subprocess.STDOUT) @@ -649,7 +667,7 @@ def run_model(self, logfile="out.log", path=None, args=None, container_path=None stopt = time.time() self.exec_time = stopt - startt if returncode.returncode == 0: - print(f"OGS finished with project file {self.project_file}.") + print(f"OGS finished with project file {self.prjfile}.") print(f"Execution took {self.exec_time} s") else: print(f"Error code: {returncode.returncode}") @@ -707,7 +725,7 @@ def write_input(self, keep_includes=False): ET.indent(self.tree, space=" ") if self.verbose is True: display.Display(self.tree) - self.tree.write(self.project_file, + self.tree.write(self.prjfile, encoding="ISO-8859-1", xml_declaration=True, pretty_print=True) @@ -716,7 +734,7 @@ def write_input(self, keep_includes=False): ET.indent(self.tree, space=" ") if self.verbose is True: display.Display(self.tree) - self.tree.write(self.project_file, + self.tree.write(self.prjfile, encoding="ISO-8859-1", xml_declaration=True, pretty_print=True)