gel.geometry

  1"""Interface to hydrogel geometry and BC information, reading required files"""
  2from .header import *
  3
  4from numbers import Number as NumberType
  5from itertools import product
  6from scipy.spatial import distance_matrix
  7
  8from .helper import *
  9
 10
 11def _give_me_disp(displacement, V, nodes):
 12    # Init
 13    u_meas = Function(V)
 14
 15    coords = V.tabulate_dof_coordinates()
 16
 17    under_const = np.zeros(len(coords))
 18    for i in range(len(coords)):
 19        this_coord = coords[i]
 20        dist_mat = distance_matrix(this_coord.reshape((1,3)), nodes)
 21        if np.min(dist_mat) < 1e-6:
 22            disp = displacement[np.argmin(dist_mat)]
 23            under_const[i] = disp[i % 3]
 24    u_meas.vector().set_local(under_const)
 25    u_meas.vector().apply("")
 26
 27    return u_meas
 28
 29
 30class _DisplacementBasedSubdomain(SubDomain):
 31
 32    _tol = 1e-14
 33
 34    def __init__(self, u, lower_bound, upper_bound):
 35        self.u = u
 36        self.lower_bound = lower_bound
 37        self.upper_bound = upper_bound
 38        super().__init__()
 39
 40    def inside(self, x, on_boundary):
 41        u_here = np.linalg.norm(self.u(x))
 42        return (
 43            (self.lower_bound - self._tol < u_here)
 44            and (u_here < self.upper_bound + self._tol)
 45        )
 46
 47
 48class Geometry:
 49
 50    CUBE = 201
 51    """Subdomain marker for far edges of gel box"""
 52    CELL = 202
 53    """Subdomain marker for inner cell surface"""
 54    UNDETECTABLE_U = 250
 55    """Subdomain marker for volume outside event horizon, when applicable"""
 56    DETECTABLE_U = 251
 57    """Subdomain marker for volume within event horizon, when applicable"""
 58
 59    def __init__(
 60            self,
 61            data_directory,
 62            xdmf_name="reference_domain",
 63            suppress_cell_bc=False,
 64            u_magnitude_subdomains_file=None,
 65            detectable_u_cutoff=0.38, # microns
 66            bci_data=None,
 67            bco_data=None
 68        ):
 69        """
 70        Object that contains data pertaining to geometry and BCs.
 71
 72        * `data_directory`: str path to a directory containing required
 73        files listed below.
 74        * `xdmf_name`: str name of an .xdmf file in `data_directory`
 75        containing FEniCS-compatible geometry with subdomain markers
 76        and the variant with boundary subdomain markers.
 77        * `suppress_cell_bc`: bool, enables removing BC on the cell
 78        surface for Neumann BC
 79        * `u_magnitude_subdomains_file`: str path to full-shape .xdmf
 80        file with displacements labeled "u" from which the event horizon
 81        is determined (None => no event horizon used)
 82        * `detectable_u_cutoff`: float cutoff value for event horizon
 83        if such functionality enabled by supplying a
 84        `u_magnitude_subdomains_file`
 85        * `bci_data`: str or None. Describes the inner cell surface
 86        boundary condition. If None, looks inside `data_directory` for
 87        "bci.vtk" with point_data "u" on cell surface. Otherwise, looks
 88        for same under path to a .vtk file.
 89        * `bco_data`: str or None. Describes the outer box boundary
 90        condition. If None, sets the BC to 0 displacement. If a str,
 91        interprets as a path to a .vtk file with displacements in "u"
 92        point_data.
 93
 94        Reads relevant data from files in `data_directory`.
 95        Required files:
 96        * {xdmf_name}.xdmf
 97        * {xdmf_name}.h5
 98        * {xdmf_name}_boundaries.xdmf
 99        * {xdmf_name}_boundaries.h5
100        * `bci_data` .vtk file
101        * If `bco_data` is not None, `bco_data` .vtk file
102        """
103        # Helper to gel relative path
104        rp = ghrp(data_directory)
105
106        # Gel Volume Mesh
107        mesh = Mesh()
108        with XDMFFile(rp(f"{xdmf_name}.xdmf")) as infile:
109            infile.read(mesh)
110
111        # Boundary info
112        mvc = MeshValueCollection("size_t", mesh, 2)
113        with XDMFFile(rp(f"{xdmf_name}_boundaries.xdmf")) as infile:
114            infile.read(mvc, "boundaries") 
115        boundaries = cpp.mesh.MeshFunctionSizet(mesh, mvc)
116
117        # Function Spaces
118        V0 = FunctionSpace(mesh, "Lagrange", 1)
119        V = VectorFunctionSpace(mesh, "Lagrange", 1)
120        VC = TensorFunctionSpace(mesh, "CG", degree=1, shape=(3,3))
121        self._DG0 = None
122
123        # Create subdomains if requested, create volumetric measure
124        if u_magnitude_subdomains_file is not None:
125            # Load in u
126            u = Function(V)  
127
128            u_file = XDMFFile(u_magnitude_subdomains_file)
129            u_file.read_checkpoint(u, "u", 0)
130
131            u.set_allow_extrapolation(True)
132            u = interpolate(u, V)
133            self.u_target = u
134
135            # Create subdomains
136            regions = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
137            
138            # Chose 40 microns since near size of cell
139            self.detectable_u_cutoff = detectable_u_cutoff
140            detectable_subdomain = _DisplacementBasedSubdomain(
141                u,
142                detectable_u_cutoff,
143                40
144            )
145
146            # Manually assign values to region due to need for alternative
147            vertex_coords = mesh.coordinates()
148            for ci, cell in enumerate(mesh.cells()):
149                point_array = vertex_coords[cell]
150
151                in_detectable_u = False
152                for point_coord in point_array:
153                    x = Point((*point_coord,))
154                    if detectable_subdomain.inside(x, False):
155                        # Only need a single point inside to count
156                        # Note default .mark behavior wants all vertices inside
157                        # and some other unknown criteria (or has a bug)
158                        in_detectable_u = True
159                        break
160
161                regions.array()[ci] = (self.DETECTABLE_U if in_detectable_u 
162                                       else self.UNDETECTABLE_U)
163
164            dx = Measure("dx", domain=mesh, subdomain_data=regions)
165        else:
166            dx = Measure("dx", domain=mesh)
167
168        # Surface measure
169        ds = Measure("ds", domain=mesh, subdomain_data=boundaries)
170
171        #
172        # Boundary Conditions
173        #
174        # Outer
175        zero = Constant((0.0, 0.0, 0.0))
176
177        if bco_data is None:
178            # Default case: 0
179            outer_bc = DirichletBC(V, zero, boundaries, 201)
180        else:
181            outer_mesh = meshio.read(bco_data)
182            outer_surface_nodes = outer_mesh.points
183            displacements = outer_mesh.point_data["u"]
184            outer_bf = _give_me_disp(displacements, V, outer_surface_nodes)
185
186            outer_bc = DirichletBC(V, outer_bf, boundaries, 201)
187
188        bcs = [outer_bc]
189
190        # Inner
191        if not suppress_cell_bc:
192            if bci_data is None:
193                # Default case: look for bci.vtk
194                bci_data = rp("bci.vtk")
195
196            surf_mesh = meshio.read(bci_data)
197            surface_nodes1 = surf_mesh.points
198            displacements = surf_mesh.point_data["u"]
199
200            self.cell_vertices = surface_nodes1
201
202            self.bf = _give_me_disp(
203                displacements,
204                V,
205                surface_nodes1
206            )
207            self.scalar = 0.0
208
209            inner_bc = DirichletBC(V, zero, boundaries, 202)
210            bcs.append(inner_bc)
211
212        # Set to internal variables
213        self.V0 = V0
214        """Scalar 1st order Lagrange function space on gel mesh"""
215        self.V = V
216        """Vector 1st order Lagrange function space on gel mesh"""
217        self.VC = VC
218        """2nd rank tensor 1st order Lagrange function space on gel mesh"""
219
220        self.dx = dx
221        """Volumetric measure on gel mesh with subdomain data"""
222        self.ds = ds
223        """Surface measure on gel mesh boundaries with subdomain data"""
224
225        self.bcs = bcs
226        """List of FEniCS boundary conditions. 0 is outer, 1 is inner
227        (if present)
228        """
229
230        self.mesh = mesh
231        """Corresponding FEniCS mesh object."""
232
233        self.boundaries = boundaries
234        """Boundary MeshFunctionSizet with surface tags"""
235        self._suppress_cell_bc = suppress_cell_bc
236
237    def output_regions(self, filename):
238        """Writes subdomain data to file under path filename."""
239        regions = self.dx.subdomain_data()
240        File(filename) << regions
241
242    def update_bcs(self):
243        """Updates `bf` according to `scalar` for load stepping.
244
245        Side effects: inner cell surface boundary condition `bcs`
246        updates according to float `scalar`
247        """
248        if not self._suppress_cell_bc:
249            new_bc = Function(self.V)
250            new_bc.vector().set_local(self.bf.vector().get_local()*self.scalar)
251            new_bc.vector().apply("")
252            self.bcs[1] = DirichletBC(self.V, new_bc, self.boundaries, 202)
253
254    @property
255    def DG0(self):
256        """Element-wise basis/function space on the gel mesh"""
257        if self._DG0 is None:
258            self._DG0 = FunctionSpace(self.mesh, "DG", 0)
259        return self._DG0