some initial paths towards model free stuff.

.gitignore

@@ -26,3 +26,5 @@ pinn*.log
 sierra_batch_script_*
 rocm_venv/
 venv/
+dd_solver.dat
+krylov_solver.dat

pancax/__init__.py

@@ -10,6 +10,7 @@
     BlatzKo, \
     Gent, \
     Hencky, \
+    InputPolyConvexPotential, \
     NeoHookean, \
     Swanson, \
     PronySeries, \
@@ -106,6 +107,7 @@
     "BlatzKo",
     "Gent",
     "Hencky",
+    "InputPolyConvexPotential",
     "NeoHookean",
     "Swanson",
     "PronySeries",

pancax/constitutive_models/__init__.py

@@ -6,6 +6,7 @@
     BlatzKo, \
     Gent, \
     Hencky, \
+    InputPolyConvexPotential, \
     NeoHookean, \
     Swanson
 from .mechanics.hyperviscoelasticity import \
@@ -26,6 +27,7 @@
     "BlatzKo",
     "Gent",
     "Hencky",
+    "InputPolyConvexPotential",
     "NeoHookean",
     "Swanson",
     # Hyperviscoelasticity

pancax/constitutive_models/mechanics/hyperelasticity/__init__.py

@@ -1,7 +1,15 @@
 from .blatz_ko import BlatzKo
 from .gent import Gent
 from .hencky import Hencky
+from .input_polyconvex_potential import InputPolyConvexPotential
 from .neohookean import NeoHookean
 from .swanson import Swanson
 
-__all__ = ["BlatzKo", "Gent", "Hencky", "NeoHookean", "Swanson"]
+__all__ = [
+    "BlatzKo",
+    "Gent",
+    "Hencky",
+    "InputPolyConvexPotential",
+    "NeoHookean",
+    "Swanson"
+]

pancax/constitutive_models/mechanics/hyperelasticity/input_polyconvex_potential.py

@@ -0,0 +1,53 @@
+from ..base import HyperelasticModel
+from ....networks.input_polyconvex_nn import InputPolyconvexNN
+from ...properties import Property
+import equinox as eqx
+import jax
+import jax.numpy as jnp
+
+
+class InputPolyConvexPotential(HyperelasticModel):
+    bulk_modulus: Property
+    network: InputPolyconvexNN
+
+    def __init__(
+        self,
+        bulk_modulus: Property,
+        key: jax.random.PRNGKey
+    ):
+        self.bulk_modulus = bulk_modulus
+        self.network = InputPolyconvexNN(
+            n_convex=2,
+            n_inputs=2,
+            n_outputs=1,
+            activation_x=jax.nn.softplus,
+            activation_y=jax.nn.softplus,
+            key=key
+        )
+
+    def energy(self, grad_u, theta, state_old, dt):
+        K = self.bulk_modulus
+        J = self.jacobian(grad_u)
+        I1_bar = self.I1_bar(grad_u)
+        I2_bar = self.I2_bar(grad_u)
+        invariants = jnp.array([I1_bar, I2_bar])
+
+        W_vol = 0.5 * K * (0.5 * (J**2 - 1) - jnp.log(J))
+        W_dev_nn = self.network(invariants)
+        W_dev_0 = self.network(jnp.array([3., 3.]))
+        # W_dev_s = 2. * jnp.sum(jax.grad(self.network))(invariants) * (J - 1.)
+        grad_W_dev_nn = jax.grad(lambda x: self.network(x)[0])(
+            jnp.array([3., 3.])
+        )
+        # TODO probably need to add dpsi/dJ to the term below
+        W_dev_s = 2. * (grad_W_dev_nn[0] + 2. * grad_W_dev_nn[1]) * (J - 1.)
+
+        return W_vol + (W_dev_nn - W_dev_0 - W_dev_s), state_old
+
+    def parameter_enforcement(self):
+        self = eqx.tree_at(
+            lambda x: x.network,
+            self,
+            self.network.parameter_enforcement()
+        )
+        return self

scripts/extract_full_field_mechanics_data.py

@@ -0,0 +1,91 @@
+from pancax.fem.read_exodus_mesh import _read_node_sets
+import argparse
+import netCDF4 as nc
+import numpy as np
+import pandas
+
+
+parser = argparse.ArgumentParser(
+    description="A script to full-field data "
+    "responses from Sierra Solid Mechanics."
+)
+
+parser.add_argument(
+    "--csv-file",
+    help="name of a file to write data to in CSV format"
+)
+parser.add_argument(
+    "--exodus-file",
+    help="The name of the exodus file"
+)
+parser.add_argument(
+    "--nodal-variables",
+    help="name of the internal force variable"
+)
+parser.add_argument(
+    "--nodeset",
+    help="name of the node set to read data from"
+)
+
+
+args = parser.parse_args()
+
+csv_file = args.csv_file
+exo_file = args.exodus_file
+nodal_var_names = args.nodal_variables.split(",")
+nset_name = args.nodeset
+
+print(f"Nodal var names = {nodal_var_names}")
+
+
+def extract_nodal_data(exo_file, nset_name, var_names):
+    vars = []
+    with nc.Dataset(exo_file, "r") as data:
+        nsets = _read_node_sets(data)
+        node_var_names = data.variables["name_nod_var"]
+        node_var_names.set_auto_mask(False)
+        node_var_names = [
+            b"".join(c).decode("UTF-8").rstrip() for c in node_var_names[:]
+        ]
+
+        assert nset_name in nsets.keys(), \
+            f"Available nodesets are {nsets.keys()}"
+        for var_name in var_names:
+            assert var_name in node_var_names
+
+        nset = nsets[nset_name]
+        times = data.variables["time_whole"][:]
+        coords = np.vstack((
+            data.variables["coordx"][nset],
+            data.variables["coordy"][nset],
+            data.variables["coordz"][nset]
+        ))
+
+        for var_name in var_names:
+            var_index = node_var_names.index(var_name)
+            var_name = f"vals_nod_var{var_index + 1}"
+            var_all_times = data.variables[var_name]
+            temp_var = []
+            for n, time in enumerate(times):
+                var_vals = np.array(var_all_times[n, nset])
+                temp_var.append(var_vals)
+            temp_var = np.hstack(temp_var)
+            vars.append(temp_var)
+
+        n_time_steps = len(times)
+        n_nodes = int(len(vars[0]) / n_time_steps)
+        coord_vals = []
+        time_vals = []
+        for n, time in enumerate(times):
+            coord_vals.append(coords)
+            time_vals.append(time * np.ones((n_nodes,)))
+        coord_vals = np.hstack(coord_vals)
+        time_vals = np.hstack(time_vals)
+
+    return np.vstack((time_vals, coord_vals, *vars)).T
+
+
+df_data = extract_nodal_data(exo_file, nset_name, nodal_var_names)
+df_columns = ["t", "x", "y", "z", *nodal_var_names]
+df = pandas.DataFrame(df_data, columns=df_columns)
+df.to_csv(csv_file, index=False)

scripts/extract_global_mechanics_data.py

@@ -0,0 +1,93 @@
+from pancax.fem.read_exodus_mesh import _read_node_sets
+import argparse
+import netCDF4 as nc
+import numpy as np
+import pandas
+
+
+parser = argparse.ArgumentParser(
+    description="A script to extract global force displacement "
+    "responses from Sierra Solid Mechanics."
+)
+
+parser.add_argument(
+    "--csv-file",
+    help="name of a file to write data to in CSV format"
+)
+parser.add_argument(
+    "--displacement-variable",
+    help="name of the displacment variable"
+)
+parser.add_argument(
+    "--exodus-file",
+    help="The name of the exodus file"
+)
+parser.add_argument(
+    "--force-variable",
+    help="name of the internal force variable"
+)
+parser.add_argument(
+    "--nodeset",
+    help="name of the node set to read data from"
+)
+
+
+args = parser.parse_args()
+
+csv_file = args.csv_file
+disp_var_name = args.displacement_variable
+exo_file = args.exodus_file
+force_var_name = args.force_variable
+nset_name = args.nodeset
+
+
+def extract_data(exo_file, nset_name, disp_var_name, force_var_name):
+    displacements = []
+    forces = []
+    with nc.Dataset(exo_file, "r") as data:
+        nsets = _read_node_sets(data)
+        node_var_names = data.variables["name_nod_var"]
+        node_var_names.set_auto_mask(False)
+        node_var_names = [
+            b"".join(c).decode("UTF-8").rstrip() for c in node_var_names[:]
+        ]
+
+        assert nset_name in nsets.keys(), \
+            f"Available nodesets are {nsets.keys()}"
+        assert disp_var_name in node_var_names
+        assert force_var_name in node_var_names
+
+        nset = nsets[nset_name]
+        disp_var_index = node_var_names.index(disp_var_name)
+        force_var_index = node_var_names.index(force_var_name)
+
+        disp_var_name = f"vals_nod_var{disp_var_index + 1}"
+        force_var_name = f"vals_nod_var{force_var_index + 1}"
+
+        times = data.variables["time_whole"][:]
+
+        disp_all_times = data.variables[disp_var_name]
+        force_all_times = data.variables[force_var_name]
+
+        for n, time in enumerate(times):
+            disp = np.array(disp_all_times[n, nset])
+            force = np.array(force_all_times[n, nset])
+
+            displacements.append(np.mean(disp))
+            forces.append(np.sum(force))
+
+    times = np.array(times)
+    displacements = np.array(displacements)
+    forces = np.array(forces)
+
+    return times, displacements, forces
+
+
+times, displacements, forces = \
+    extract_data(exo_file, nset_name, disp_var_name, force_var_name)
+
+# save sierra exo file to csv file
+df_data = np.array([times, displacements, forces]).T
+df_columns = ["times", "disps", "forces"]
+df = pandas.DataFrame(df_data, columns=df_columns)
+df.to_csv(csv_file, index=False)

examples/forward_problems/mechanics/plot_force_disp.py → scripts/plot_force_disp.py

@@ -5,10 +5,12 @@
 import numpy as np
 import pandas
 
-parser = argparse.ArgumentParser(description="A script to plot force displacement curves.")
+parser = argparse.ArgumentParser(
+    description="A script to plot force displacement curves."
+)
 
 parser.add_argument(
-    "--displacement-variable", 
+    "--displacement-variable",
     help="name of the displacment variable"
 )
 parser.add_argument(
@@ -24,7 +26,7 @@
     help="name of the Sierra node set to read data from"
 )
 parser.add_argument(
-    "--sierra-exodus-file", 
+    "--sierra-exodus-file",
     help="The name of the Sierra exodus file."
 )
 parser.add_argument(
@@ -49,9 +51,12 @@ def extract_data(exo_file, nset_name, disp_var_name, force_var_name):
         nsets = _read_node_sets(data)
         node_var_names = data.variables["name_nod_var"]
         node_var_names.set_auto_mask(False)
-        node_var_names = [b"".join(c).decode("UTF-8").rstrip() for c in node_var_names[:]]
+        node_var_names = [
+            b"".join(c).decode("UTF-8").rstrip() for c in node_var_names[:]
+        ]
 
-        assert nset_name in nsets.keys(), f"Available nodesets are {nsets.keys()}"
+        assert nset_name in nsets.keys(), \
+            f"Available nodesets are {nsets.keys()}"
         assert disp_var_name in node_var_names
         assert force_var_name in node_var_names
 
@@ -80,8 +85,13 @@ def extract_data(exo_file, nset_name, disp_var_name, force_var_name):
 
     return times, displacements, forces
 
-p_times, p_displacements, p_forces = extract_data(pinn_exo_file, pinn_nset_name, disp_var_name, force_var_name)
-s_times, s_displacements, s_forces = extract_data(sierra_exo_file, sierra_nset_name, disp_var_name, force_var_name)
+
+p_times, p_displacements, p_forces = \
+    extract_data(pinn_exo_file, pinn_nset_name, disp_var_name, force_var_name)
+s_times, s_displacements, s_forces = \
+    extract_data(
+        sierra_exo_file, sierra_nset_name, disp_var_name, force_var_name
+    )
 
 # save sierra exo file to csv file
 df_data = np.array([p_times, p_displacements, p_forces]).T

scripts/run_linter.sh

@@ -1,2 +1,5 @@
 flake8 pancax
 flake8 test
+flake8 scripts/extract_full_field_mechanics_data.py
+flake8 scripts/extract_global_mechanics_data.py
+flake8 scripts/plot_force_disp.py