Add multiple laws, parameter method

python/comfe/cdm.py

@@ -100,12 +100,13 @@ def __init__(
         f_ext: Callable | None,
         bcs: list[df.fem.DirichletBCMetaClass],
         M: df.fem.Function,
-        rust_model: RustConstitutiveModel,
+        rust_model: RustConstitutiveModel | list[RustConstitutiveModel],
         quadrature_rule: QuadratureRule,
         additional_output: list[str] | None = None,
         calculate_total_energy: bool = False,
         calculate_bulk_viscosity: bool = False,
         monitor_distortion: bool = False,
+        cells: list[np.ndarray[np.int32]] | None = None,
     ) -> None:
         # self.del_t = None
         v = df.fem.Function(function_space, name="Velocity")
@@ -122,11 +123,18 @@ def __init__(
         else:
             distortion_expr = None
 
+        ips = [] if cells is not None else None
+        cells = [] if cells is None else cells
+        for cell_array in cells:
+            number_of_ips = quadrature_rule.weights.size
+            ips_array = [cell_array * number_of_ips + i for i in range(number_of_ips)]
+            ips.append(np.concatenate(ips_array, dtype=np.uint64, axis=None))
+
         model = ConstitutiveModel(
             rust_model,
             quadrature_rule,
             function_space.mesh,
-            None,
+            ips,
             additional_output,
         )
 
@@ -654,7 +662,7 @@ def __init__(
         t0: float,
         f_ext: Callable | None,
         bcs: list[df.fem.DirichletBCMetaClass],
-        rust_model: RustConstitutiveModel,
+        rust_model: RustConstitutiveModel | list[RustConstitutiveModel],
         parameters: dict[str, float],
         nonlocal_parameters: dict[str, float],
         Q_local: str,
@@ -669,12 +677,23 @@ def __init__(
         mass_nonlocal: df.fem.Function | None = None,
         calculate_bulk_viscosity: bool = False,
         nonlocal_initial_config: bool = True,
+        cells: list[np.ndarray[np.int32]] | None = None,
     ) -> None:
-        mass_mechanics = (
-            mass_mechanics
-            if mass_mechanics is not None
-            else diagonal_mass(velocity_space, parameters["rho"], invert=True)
-        )
+        if mass_mechanics is not None:
+            rho_space = df.fem.FunctionSpace(velocity_space.mesh, ("DG", 0))
+            rho = df.fem.Function(rho_space)
+            cells = cells if cells is not None else [np.arange(rho.x.array.size, dtype=np.int32)]
+            rho_list = (
+                [model.parameters()["RHO"] for model in rust_model] if cells is not None else [parameters["rho"]]
+            )
+            assert len(cells) == len(rho_list)
+
+            for cells_i, rho_i in zip(cells, rho_list):
+                rho.x.array[cells_i] = rho_i
+                rho.x.scatter_forward()
+
+            mass_mechanics = diagonal_mass(velocity_space, rho, invert=True)
+
         mass_nonlocal = (
             mass_nonlocal
             if mass_nonlocal is not None or nonlocal_solver == ImplicitNonlocal

python/comfe/laws.py

@@ -1,4 +1,5 @@
 from abc import ABC, abstractmethod
+from functools import reduce
 
 import dolfinx as df
 import numpy as np
@@ -35,30 +36,39 @@ def update(self) -> None:
 
 
 class ConstitutiveModel(BaseModel):
-    rs_object: RustConstitutiveModel
+    rs_object: list[RustConstitutiveModel]
     input: dict[str, df.fem.Function]
     output: dict[str, df.fem.Function]
-    ips: np.ndarray[np.uint64] | None = None
+    ips: list[np.ndarray[np.uint64]] | None = None
     spaces: dict[int | tuple[int, int], df.fem.FunctionSpace]
 
     class Config:
         arbitrary_types_allowed = True
 
     def __init__(
         self,
-        model: RustConstitutiveModel,
+        model: RustConstitutiveModel | list[RustConstitutiveModel],
         rule: QuadratureRule,
         mesh: df.mesh.Mesh,
-        ips: np.ndarray[np.uint64] | None = None,
+        ips: np.ndarray[np.uint64] | list[np.ndarray[np.uint64]] | None = None,
         additional_variables: list[str] | None = None,
     ) -> None:
-        input, output, spaces = ceate_input_and_output(model, rule, mesh, None, additional_variables)
+        if isinstance(model, RustConstitutiveModel):
+            model = [model]
+
+        # checks that all models are of the same type
+        assert reduce(lambda x,y: y if type(x) is type(y) else None, model) is not None
+        input, output, spaces = ceate_input_and_output(model[0], rule, mesh, None, additional_variables)
         super().__init__(rs_object=model, input=input, output=output, ips=ips, spaces=spaces)
 
     def evaluate(self, del_t=1.0) -> None:
         input = {key: value.vector.array for key, value in self.input.items()}
         output = {key: value.vector.array for key, value in self.output.items()}
-        self.rs_object.evaluate(del_t, input, output)
+        if len(self.rs_object) > 1:
+            for model, ips in zip(self.rs_object, self.ips):
+                model.evaluate_some(del_t, ips, input, output)
+        else:
+            self.rs_object[0].evaluate(del_t, input, output)
 
     def evaluate_some(self, del_t=1.0) -> None:
         input = {key: value.vector.array for key, value in self.input.items()}

src/gradient_jh2.rs

@@ -315,4 +315,35 @@ impl ConstitutiveModel for GradientJH23D {
             (Q::InternalHeatingEnergy, QDim::Scalar),
         ])
     }
+
+    fn parameters(&self)->HashMap<String, f64> {
+        HashMap::from(
+            [
+                ("RHO".to_string(), self.parameters.RHO),
+                ("SHEAR_MODULUS".to_string(), self.parameters.SHEAR_MODULUS),
+                ("A".to_string(), self.parameters.A),
+                ("B".to_string(), self.parameters.B),
+                ("C".to_string(), self.parameters.C),
+                ("M".to_string(), self.parameters.M),
+                ("N".to_string(), self.parameters.N),
+                ("EPS0".to_string(), self.parameters.EPS0),
+                ("T".to_string(), self.parameters.T),
+                ("SIGMAHEL".to_string(), self.parameters.SIGMAHEL),
+                ("PHEL".to_string(), self.parameters.PHEL),
+                ("D1".to_string(), self.parameters.D1),
+                ("D2".to_string(), self.parameters.D2),
+                ("K1".to_string(), self.parameters.K1),
+                ("K2".to_string(), self.parameters.K2),
+                ("K3".to_string(), self.parameters.K3),
+                ("BETA".to_string(), self.parameters.BETA),
+                ("EFMIN".to_string(), self.parameters.EFMIN),
+                ("DMAX".to_string(), self.parameters.DMAX),
+                ("E_F".to_string(), self.parameters.E_F),
+                ("E_0".to_string(), self.parameters.E_0),
+                ("LOCAL_SOUND_SPEED".to_string(), self.parameters.LOCAL_SOUND_SPEED),
+                ("HARDENING".to_string(), self.parameters.HARDENING),
+                ("M_OVERNONLOCAL".to_string(), self.nonlocal_parameters.m_overnonlocal),
+            ]
+        )
+    }
 }

src/hypoelasticity.rs

@@ -68,4 +68,11 @@ impl ConstitutiveModel for Hypoelasticity3D {
     fn define_optional_history(&self) -> HashMap<Q, QDim> {
         HashMap::from([])
     }
+
+    fn parameters(&self)->HashMap<String, f64> {
+        HashMap::from([
+            ("mu".to_string(), self.mu),
+            ("lambda".to_string(), self.lambda),
+        ])
+    }
 }

src/interfaces.rs

@@ -311,6 +311,7 @@ pub trait ConstitutiveModel {
             return Err("There are inconsistencies in input and output sizes.");
         }
     }
+    fn parameters(&self)->HashMap<String, f64>;
 }
 
 

src/jh2.rs

@@ -329,4 +329,32 @@ impl ConstitutiveModel for JH23D {
             (Q::InternalHeatingEnergy, QDim::Scalar),
         ])
     }
+
+    fn parameters(&self)->HashMap<String, f64> {
+        HashMap::from([
+            ("RHO".to_string(), self.parameters.RHO),
+            ("SHEAR_MODULUS".to_string(), self.parameters.SHEAR_MODULUS),
+            ("A".to_string(), self.parameters.A),
+            ("B".to_string(), self.parameters.B),
+            ("C".to_string(), self.parameters.C),
+            ("M".to_string(), self.parameters.M),
+            ("N".to_string(), self.parameters.N),
+            ("EPS0".to_string(), self.parameters.EPS0),
+            ("T".to_string(), self.parameters.T),
+            ("SIGMAHEL".to_string(), self.parameters.SIGMAHEL),
+            ("PHEL".to_string(), self.parameters.PHEL),
+            ("D1".to_string(), self.parameters.D1),
+            ("D2".to_string(), self.parameters.D2),
+            ("K1".to_string(), self.parameters.K1),
+            ("K2".to_string(), self.parameters.K2),
+            ("K3".to_string(), self.parameters.K3),
+            ("BETA".to_string(), self.parameters.BETA),
+            ("EFMIN".to_string(), self.parameters.EFMIN),
+            ("DMAX".to_string(), self.parameters.DMAX),
+            ("E_F".to_string(), self.parameters.E_F),
+            ("E_0".to_string(), self.parameters.E_0),
+            ("LOCAL_SOUND_SPEED".to_string(), self.parameters.LOCAL_SOUND_SPEED),
+            ("HARDENING".to_string(), self.parameters.HARDENING),
+        ])
+    }
 }

src/lib.rs

@@ -3,7 +3,7 @@ use std::collections::HashMap;
 
 use crate::interfaces::{ConstitutiveModel, QDim, QValueInput, QValueOutput, Q};
 use crate::jh2::JH23D;
-use crate::jhr::JHR3D;
+//use crate::jhr::JHR3D;
 //use crate::jh_concrete::JHConcrete3D;
 use crate::generic_jh2::GenericJH23D;
 use crate::gradient_jh2::GradientJH23D;
@@ -22,7 +22,7 @@ pub mod jh2;
 //pub mod jh_concrete;
 pub mod generic_jh2;
 pub mod gradient_jh2;
-pub mod jhr;
+//pub mod jhr;
 pub mod hypoelasticity;
 pub mod smallstrain;
 pub mod stress_strain;
@@ -237,6 +237,14 @@ macro_rules! impl_constitutive_model {
                 }
                 Ok(history_py.into())
             }
+            fn parameters(&self, py: Python) -> PyResult<PyObject> {
+                let parameters_py = PyDict::new(py);
+                let parameters_rs = self.model.parameters();
+                for (key, value) in parameters_rs.iter() {
+                    parameters_py.set_item(key, value)?;
+                }
+                Ok(parameters_py.into())
+            }
         }
         $m.add_class::<$name>()?;
     };

src/smallstrain/linear_elastic.rs

@@ -5,6 +5,8 @@ use std::collections::HashMap;
 #[derive(Debug)]
 pub struct LinearElastic3D {
     pub D: SMatrix<f64, 6, 6>,
+    pub E: f64,
+    pub nu: f64,
 }
 
 
@@ -53,7 +55,7 @@ impl ConstitutiveModel for LinearElastic3D {
             0.0,
             2.0 * mu,
         );
-        Some(Self { D: D })
+        Some(Self { D: D , E: *E, nu: *nu})
     }
     fn evaluate_ip(&self, ip: usize, _del_t: f64, input: &QValueInput, output: &mut QValueOutput) {
         let strain = input.get_vector::<{Q::MandelStrain.size()}>(Q::MandelStrain, ip);
@@ -74,5 +76,11 @@ impl ConstitutiveModel for LinearElastic3D {
             (Q::MandelTangent, Q::MandelTangent.q_dim()),
         ])
     }
+    fn parameters(&self)-> HashMap<String, f64>{
+        HashMap::from([
+            ("E".to_string(), self.E),
+            ("nu".to_string(), self.nu),
+        ])
+    }
 
 }