hermit.fenics.ops ================= .. py:module:: hermit.fenics.ops .. autoapi-nested-parse:: The FEniCSx custom operations. Called only from ``hermit._solve`` / ``hermit.outputs`` / ``hermit.transfer``. ``ShellSolveOp`` -- implicit; assemble R/J, linear (or Newton) solve, adjoint. ``ShellScalarFormsOp`` (scalar output forms) -- explicit, dense Jacobians. ``ShellFieldFormsOp`` (field outputs -- one UFL expression per :class:`FieldSpec` into its target space, by ``project`` / ``interpolate`` / ``average`` / ``midpoint``) -- explicit, reverse-mode only (one cotangent solve or matvec per field, never a dense Jacobian). ``project`` / ``average`` / ``midpoint`` keep mesh-coordinate derivatives; ``interpolate`` does not. Port of ``StateOperation`` / ``OutputOperation`` / ``OutputFieldOperation`` (femo dev_coupling), self-contained (no ``FEA``). The shell tangent is symmetric, so one MUMPS factorization of ``dR/dw`` serves both the forward solve and the adjoint. .. !! processed by numpydoc !! Classes ------- .. autoapisummary:: hermit.fenics.ops.FieldSpec hermit.fenics.ops.ShellFieldFormsOp hermit.fenics.ops.ShellScalarFormsOp hermit.fenics.ops.ShellSolveOp Functions --------- .. autoapisummary:: hermit.fenics.ops.orientation_cos_sin_spec hermit.fenics.ops.strain_fields Module Contents --------------- .. py:class:: FieldSpec One field output: a UFL expression + how to represent it as FE coefficients. ``space`` is ``(family, degree)`` (blocked by ``n_components`` when the op builds the FunctionSpace). ``method``: * ``"project"`` (default) -- L2 projection ``M c = b``. Any space; one mass factorization; full derivatives including mesh coordinates. * ``"interpolate"`` -- collocate ``expr`` at the space's nodal points. Any Lagrange space; cheaper (no solve); exact where ``expr`` already lies in the space. Has **no** mesh-coordinate (shape) derivative -- use it only when the geometry is fixed. * ``"average"`` / ``"midpoint"`` -- ``(1/vol_K) * integral_K(expr)``, i.e. L2 projection onto ``("DG", 0)`` (diagonal mass -> a divide, no solve); ``"average"`` uses full quadrature, ``"midpoint"`` a single centroid point. DG0 only; full derivatives including mesh coordinates. .. !! processed by numpydoc !! .. py:attribute:: args :type: tuple :value: ('disp_solid',) .. py:attribute:: expr :type: object .. py:attribute:: frame :type: str :value: 'local' .. py:attribute:: kind :type: str :value: 'scalar' .. py:attribute:: method :type: str :value: 'project' .. py:attribute:: n_components :type: int .. py:attribute:: space :type: tuple :value: ('DG', 2) .. py:class:: ShellFieldFormsOp(pde, fields: dict, differentiable_geometry: bool = False, coefficients=None, quadrature_degree=4) Bases: :py:obj:`csdl_alpha.CustomExplicitOperation` Field outputs -- several UFL expressions -> FE coefficients, in one op. ``fields`` maps output name -> :class:`FieldSpec`. Each output is ``(n_scalar_dofs, n_components)`` in that field's target space, FE dof order (wrap in :class:`hermit._field.Field` for user ordering, point evaluation and frame changes). Reverse-mode only (``compute_jacvec_product``), like femo's ``OutputFieldOperation``: one cotangent solve / matvec per field, never a dense Jacobian. ``method="interpolate"`` has no ``mesh_nodes`` derivative -- the op raises if ``differentiable_geometry`` is set with an interpolated field. .. !! processed by numpydoc !! .. py:method:: compute(input_vals, output_vals) .. py:method:: compute_jacvec_product(inputs, outputs, d_inputs, d_outputs, mode) .. py:method:: evaluate(fe) .. py:attribute:: arg_names .. py:attribute:: fields .. py:attribute:: gdim .. py:attribute:: nel .. py:attribute:: pde .. py:attribute:: tdim .. py:class:: ShellScalarFormsOp(pde, forms: dict, differentiable_geometry: bool = False, orientation=None, coefficients=None) Bases: :py:obj:`csdl_alpha.CustomExplicitOperation` Assemble several scalar UFL output forms + their partials in one op. Port of femo's ``OutputOperation`` (multi-form). ``forms`` maps output name -> (ufl_form, tuple of arg names). Function args are written onto the persistent ``ShellPDE`` functions; ``mesh_nodes`` writes the mesh geometry (scratch buffer, restored after each call) and, when differentiable, contributes a ``d(form)/d(SpatialCoordinate)`` partial scattered back to user node ordering. ``orientation``, given, is ``(name, space)`` (``name`` in ``("fiber_angle", "fiber_direction")``) -- needed whenever one of ``forms``' arg-name tuples includes it (``elastic_energy`` depends on theta for an oriented composite; a strain-field form does not). Unlike the fixed-space ``_ARG_TO_PDE`` args, an orientation coefficient's space is caller-chosen, so it is resolved through ``pde.coefficient(name, space)`` (see that method) rather than a fixed ``pde`` attribute -- the same seam ``ShellSolveOp`` uses, so both ops reference the exact same persistent ``Function`` the form the caller built was differentiated against. ``coefficients``, given, is a dict ``{name: space}`` generalizing that same seam to any other arg name a caller-built form needs on a caller-chosen space (Phase F's per-load-term / arbitrary-``constitutive_space`` forms -- e.g. ``hm.solve``'s ``ShellPDE.compliance_form(loads=...)``, whose args are per-term names like ``"traction_0"`` with no fixed ``pde`` attribute at all). ``None`` (default, every pre-Phase-F call site) leaves every non-orientation arg on the fixed ``_ARG_TO_PDE`` route, unchanged. .. !! processed by numpydoc !! .. py:method:: compute(input_vals, output_vals) .. py:method:: compute_derivatives(input_vals, output_vals, derivatives) .. py:method:: evaluate(fe) .. py:attribute:: arg_names .. py:attribute:: forms .. py:attribute:: gdim .. py:attribute:: pde .. py:class:: ShellSolveOp(pde, bc, arg_names, options, form_cache=None, orientation=None, material=None, loads=None) Bases: :py:obj:`csdl_alpha.experimental.CustomImplicitOperation` ``orientation``, given, is ``(name, space)`` with ``name`` in ``("fiber_angle", "fiber_direction")`` -- the in-form ``Teps(theta)`` rotation (see ``hermit.fenics.elastic_model.ElasticModel``). ``name`` must also appear in ``arg_names`` for it to be a declared (differentiable) CSDL input. ``material``, given, is a dict ``{"A": space, "B": space, "D": space, "As": space}`` (all four): each coefficient is routed through ``pde.coefficient(name, space)`` instead of the PDE's fixed ``VABD``/``VAs`` space, so an arbitrary caller-chosen ``constitutive_space`` reaches the solve. ``None`` (default) keeps the fixed-space Functions. ``loads``, given, is a list of ``(argname, kind, space)`` triples: ``kind`` in ``"traction"`` / ``"moment"`` / ``"pressure"``, each routed through ``pde.coefficient(argname, space)`` and contributing its own residual/compliance term (see ``ElasticModel.weak_residual`` / ``ShellPDE.compliance_form``'s ``load_terms``), rather than being reduced onto the fixed ``f``/``m`` VF coefficients. Every ``argname`` must also appear in ``arg_names``. ``None`` (default) keeps the fixed-VF ``f``/``m`` mechanism; an empty list selects the per-term path with zero distributed terms (e.g. a point-load-only ``Loads``). The residual is not structurally fixed once orientation, material spaces, load composition, and multi-region penalty terms can all vary -- ``form_cache`` is keyed by a structural signature of all four, not by a single fixed string, so two structurally different compositions on the same ``ShellDomain`` never share a (wrong) compiled residual. .. !! processed by numpydoc !! .. py:method:: apply_inverse_jacobian(inputs, outputs, d_outputs, d_residuals, mode) .. py:method:: compute_jacvec_product(inputs, outputs, d_inputs, d_outputs, d_residuals, mode) .. py:method:: evaluate(fe) .. py:method:: solve_residual_equations(inputs, outputs) .. py:attribute:: Vc .. py:attribute:: arg_names .. py:attribute:: bc .. py:attribute:: dR_dw .. py:attribute:: gdim .. py:attribute:: ndof .. py:attribute:: options .. py:attribute:: orientation :value: None .. py:attribute:: pde .. py:attribute:: residual .. py:attribute:: state_name :value: 'disp_solid' .. py:function:: orientation_cos_sin_spec(pde, orientation_name, coeff_space, target_space=('DG', 0), method='average') ``FieldSpec`` for ``(cos theta, sin theta)`` of the fibre orientation angle. The angle between a global ``fiber_direction`` and each element's in-plane frame is a function of the mesh: moving nodes rotates the frame. Building the cos/sin in UFL -- the same expression the shell form itself uses -- keeps that dependence in the differentiated path, so the geometry derivative reaches the Tsai-Wu outputs. Evaluating the angle in numpy instead freezes it at the reference geometry and silently drops a term from the shape gradient. ``fiber_angle`` is element-relative and carries no geometry dependence, but it is accepted here too so both orientation kinds take one path. .. !! processed by numpydoc !! .. py:function:: strain_fields(pde, space=None, method=None, frame=None) ``name -> FieldSpec`` for the membrane strain / bending curvature / transverse shear (Voigt / engineering convention). Natural space ``("DG", 2)``, natural method ``"project"``; override with ``space`` / ``method``. ``frame`` selects the components: * ``"local"`` -- the element in-plane frame: strain / curvature as a length-3 in-plane Voigt vector, shear as a length-2 vector. Only well defined on a **discontinuous** space (adjacent elements have different in-plane frames, so a continuous space would average frame-inconsistent components at shared nodes). * ``"global"`` -- global Cartesian: strain / curvature as a length-6 symmetric tensor ``[xx, yy, zz, 2yz, 2xz, 2xy]`` (engineering shear), shear as a length-3 vector. Frame-consistent -> works on CG spaces. ``frame=None`` picks ``"local"`` for a DG space and ``"global"`` for a CG space. .. !! processed by numpydoc !!