79 add random field feature

.vscode/settings.json

@@ -0,0 +1,8 @@
+{
+    "python.testing.pytestArgs": [
+        "tests"
+    ],
+    "python.testing.unittestEnabled": false,
+    "python.testing.pytestEnabled": true,
+    "python.analysis.typeCheckingMode": "basic"
+}

environment.yml

@@ -8,6 +8,8 @@ dependencies:
   - matplotlib
   - python-gmsh
   - meshio
+  - jsonschema
+  - scipy
   # tests
   - pytest
   - coverage

src/fenicsxconcrete/finite_element_problem/base_material.py

@@ -1,11 +1,15 @@
+import importlib
+import json
 from abc import ABC, abstractmethod
 from pathlib import Path, PosixPath
 
+import jsonschema
 import pint
 
 from fenicsxconcrete.experimental_setup.base_experiment import Experiment
 from fenicsxconcrete.helper import LogMixin, Parameters
 from fenicsxconcrete.sensor_definition.base_sensor import BaseSensor
+from fenicsxconcrete.sensor_definition.sensor_schema import generate_sensor_schema
 from fenicsxconcrete.unit_registry import ureg
 
 
@@ -99,6 +103,60 @@ def delete_sensor(self) -> None:
         del self.sensors
         self.sensors = self.SensorDict()
 
+    def export_sensors_metadata(self, path: Path) -> None:
+        """Exports sensor metadata to JSON file according to the appropriate schema.
+
+        Args:
+            path : Path
+                Path where the metadata should be stored
+
+        """
+
+        sensors_metadata_dict = {"sensors": []}
+
+        for key, value in self.sensors.items():
+            sensors_metadata_dict["sensors"].append(value.report_metadata())
+            # sensors_metadata_dict[key]["name"] = key
+
+        with open(path, "w") as f:
+            json.dump(sensors_metadata_dict, f)
+
+    def import_sensors_from_metadata(self, path: Path) -> None:
+        """Import sensor metadata to JSON file and validate with the appropriate schema.
+
+        Args:
+            path : Path
+                Path where the metadata file is
+
+        """
+
+        # Load and validate
+        sensors_metadata_dict = {}
+        with open(path, "r") as f:
+            sensors_metadata_dict = json.load(f)
+        schema = generate_sensor_schema()
+        jsonschema.validate(instance=sensors_metadata_dict, schema=schema)
+
+        for sensor in sensors_metadata_dict["sensors"]:
+            # Dynamically import the module containing the class
+            module_name = "fenicsxconcrete.sensor_definition." + sensor["sensor_file"].lower()
+            module = importlib.import_module(module_name)
+
+            # Create a dictionary of keyword arguments from the remaining properties in the dictionary
+            kwargs = {
+                k: v for k, v in sensor.items() if k not in ["id", "type", "sensor_file", "units", "dimensionality"]
+            }
+
+            # Dynamically retrieve the class by its name
+            class_name = sensor["type"]
+            MySensorClass = getattr(module, class_name)
+
+            # Instantiate an object of the class with the given properties
+            sensor_i = MySensorClass(name=sensor["id"], **kwargs)
+            sensor_i.set_units(units=sensor["units"])
+
+            self.add_sensor(sensor_i)
+
     class SensorDict(dict):
         """
         Dict that also allows to access the parameter p["parameter"] via the matching attribute p.parameter

src/fenicsxconcrete/finite_element_problem/linear_elasticity_grf.py

@@ -0,0 +1,143 @@
+from pathlib import Path, PosixPath
+
+import dolfinx as df
+import numpy as np
+import pint
+import ufl
+from petsc4py.PETSc import ScalarType
+
+from fenicsxconcrete.experimental_setup.base_experiment import Experiment
+from fenicsxconcrete.finite_element_problem.linear_elasticity import LinearElasticity
+from fenicsxconcrete.gaussian_random_field import Randomfield
+
+
+class LinearElasticityGRF(LinearElasticity):
+    """Material definition for linear elasticity"""
+
+    def __init__(
+        self,
+        experiment: Experiment,
+        parameters: dict[str, pint.Quantity],
+        pv_name: str = "pv_output_full",
+        pv_path: str = None,
+    ) -> None:
+        super().__init__(experiment, parameters, pv_name, pv_path)
+        # TODO There should be more elegant ways of doing this
+        self.pv_egrf_file = Path(pv_path) / (pv_name + "egrf.xdmf")
+        self.pv_nugrf_file = Path(pv_path) / (pv_name + "nugrf.xdmf")
+
+    def setup(self) -> None:
+        self.field_function_space = df.fem.FunctionSpace(self.experiment.mesh, ("CG", 1))
+        self.lambda_ = df.fem.Function(self.field_function_space)
+        self.mu = df.fem.Function(self.field_function_space)
+
+        lame1, lame2 = self.get_lames_constants()
+        self.lambda_.vector[:] = lame1
+        self.mu.vector[
+            :
+        ] = lame2  # make this vector as a fenics constant array. Update the lame1 and lame2 in each iteration.
+
+        # define function space ets.
+        self.V = df.fem.VectorFunctionSpace(self.mesh, ("Lagrange", self.p["degree"]))  # 2 for quadratic elements
+        self.V_scalar = df.fem.FunctionSpace(self.mesh, ("Lagrange", self.p["degree"]))
+
+        # Define variational problem
+        self.u_trial = ufl.TrialFunction(self.V)
+        self.v = ufl.TestFunction(self.V)
+
+        # initialize L field, not sure if this is the best way...
+        zero_field = df.fem.Constant(self.mesh, ScalarType(np.zeros(self.p["dim"])))
+        self.L = ufl.dot(zero_field, self.v) * ufl.dx
+
+        # apply external loads
+        external_force = self.experiment.create_force_boundary(self.v)
+        if external_force:
+            self.L = self.L + external_force
+
+        body_force = self.experiment.create_body_force(self.v)
+        if body_force:
+            self.L = self.L + body_force
+
+        # boundary conditions only after function space
+        bcs = self.experiment.create_displacement_boundary(self.V)
+
+        self.a = ufl.inner(self.sigma(self.u_trial), self.epsilon(self.v)) * ufl.dx
+        self.weak_form_problem = df.fem.petsc.LinearProblem(
+            self.a,
+            self.L,
+            bcs=bcs,
+            petsc_options={"ksp_type": "preonly", "pc_type": "lu"},
+        )
+
+    # Random E and nu fields.
+    def random_field_generator(
+        self,
+        field_function_space,
+        cov_name,
+        mean,
+        correlation_length1,
+        correlation_length2,
+        variance,
+        no_eigen_values,
+        ktol,
+    ):
+        random_field = Randomfield(
+            field_function_space,
+            cov_name,
+            mean,
+            correlation_length1,
+            correlation_length2,
+            variance,
+            no_eigen_values,
+            ktol,
+        )
+        # random_field.solve_covariance_EVP()
+        return random_field
+
+    def parameters_conversion(self, lognormal_mean, lognormal_sigma):
+        from math import sqrt
+
+        normal_mean = np.log(lognormal_mean / sqrt(1 + (lognormal_sigma / lognormal_mean) ** 2))
+        normal_sigma = np.log(1 + (lognormal_sigma / lognormal_mean) ** 2)
+        return normal_mean, normal_sigma
+
+    def get_lames_constants(
+        self,
+    ):
+        # Random E and nu fields.
+        E_mean, E_variance = self.parameters_conversion(self.p["E"], 100e9)  # 3
+        Nu_mean, Nu_variance = self.parameters_conversion(self.p["nu"], 0.3)  # 0.03
+
+        self.E_randomfield = self.random_field_generator(
+            self.field_function_space, "squared_exp", E_mean, 0.3, 0.05, E_variance, 3, 0.01
+        )
+        self.E_randomfield.create_random_field(_type="random", _dist="LN")
+
+        self.nu_randomfield = self.random_field_generator(
+            self.field_function_space, "squared_exp", Nu_mean, 0.3, 0.05, Nu_variance, 3, 0.01
+        )
+        self.nu_randomfield.create_random_field(_type="random", _dist="LN")
+
+        lame1 = (self.E_randomfield.field.vector[:] * self.nu_randomfield.field.vector[:]) / (
+            (1 + self.nu_randomfield.field.vector[:]) * (1 - 2 * self.nu_randomfield.field.vector[:])
+        )
+        lame2 = self.E_randomfield.field.vector[:] / (2 * (1 + self.nu_randomfield.field.vector[:]))
+        return lame1, lame2
+
+    # TODO move this to sensor definition!?!?!
+    def pv_plot(self, t: int = 0) -> None:
+        # TODO add possibility for multiple time steps???
+        # Displacement Plot
+
+        # "Displacement.xdmf"
+        # pv_output_file
+        # TODO Look into how to unify in one file
+        with df.io.XDMFFile(self.mesh.comm, self.pv_output_file, "w") as xdmf:
+            xdmf.write_mesh(self.mesh)
+            xdmf.write_function(self.displacement)
+        with df.io.XDMFFile(self.mesh.comm, self.pv_egrf_file, "w") as xdmf:
+            xdmf.write_mesh(self.mesh)
+            xdmf.write_function(self.E_randomfield.field)
+        with df.io.XDMFFile(self.mesh.comm, self.pv_nugrf_file, "w") as xdmf:
+            xdmf.write_mesh(self.mesh)
+            xdmf.write_function(self.nu_randomfield.field)