Add Dirichlet Boundary Projectors

psydac/api/feec.py

@@ -1,3 +1,10 @@
+import numpy as np
+
+from scipy.sparse                               import dia_matrix
+
+from sympde.expr                                import integral, BilinearForm
+from sympde.topology                            import elements_of, Line, Derham
+
 from psydac.api.basic                           import BasicDiscrete
 
 from psydac.feec.derivatives                    import Derivative1D, Gradient2D, Gradient3D
@@ -26,7 +33,12 @@
 
 from psydac.fem.basic                           import FemSpace, FemLinearOperator
 from psydac.fem.vector                          import VectorFemSpace
-from psydac.linalg.basic                        import IdentityOperator
+
+from psydac.linalg.basic                        import LinearOperator, IdentityOperator
+from psydac.linalg.block                        import BlockLinearOperator
+from psydac.linalg.direct_solvers               import BandedSolver
+from psydac.linalg.kron                         import KroneckerLinearSolver, KroneckerStencilMatrix
+from psydac.linalg.stencil                      import StencilVectorSpace
 
 __all__ = ('DiscreteDeRham', 'DiscreteDeRhamMultipatch',)
 
@@ -288,6 +300,39 @@ def derivatives(self, kind='femlinop'):
         elif kind == 'linop': 
             return tuple(b_diff.linop for b_diff in self._derivatives)
 
+    #--------------------------------------------------------------------------
+    def dirichlet_projectors(self, kind='femlinop'):
+        """
+        Returns operators that apply the correct Dirichlet boundary conditions.
+
+        Parameters
+        ----------
+        kind : str
+            The kind of the projector, can be 'femlinop' or 'linop'.
+            - 'femlinop' returns a psydac FemLinearOperator (default)
+            - 'linop' returns a psydac LinearOperator
+
+        Returns
+        -------
+        d_projectors : list
+            List of <psydac.fem.basic.FemLinearOperator> or <psydac.linalg.basic.LinearOperator>
+            The Dirichlet boundary projectors of each space and in desired form.
+
+        Notes
+        -----
+        See examples/vector_potential_3d.py for a use case of these operators in LinearOperator form.
+
+        """
+        assert kind in ('femlinop', 'linop')
+
+        from psydac.fem.projectors import DirichletBoundaryProjector
+        d_projectors = [DirichletBoundaryProjector(Vh) for Vh in self.spaces[:-1]]
+
+        if kind == 'femlinop':
+            d_projectors = [FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces[:-1], d_projectors)]
+
+        return d_projectors
+
     #--------------------------------------------------------------------------
     def conforming_projectors(self, kind='femlinop', mom_pres=False, p_moments=-1, hom_bc=False):
         """
@@ -606,3 +651,36 @@ def projectors(self, *, kind='global', nquads=None):
 
         elif self.dim == 3:
             raise NotImplementedError("3D projectors are not available")
+
+    #--------------------------------------------------------------------------
+    def dirichlet_projectors(self, kind='femlinop'):
+        """
+        Returns operators that apply the correct Dirichlet boundary conditions.
+
+        Parameters
+        ----------
+        kind : str
+            The kind of the projector, can be 'femlinop' or 'linop'.
+            - 'femlinop' returns a psydac FemLinearOperator (default)
+            - 'linop' returns a psydac LinearOperator
+
+        Returns
+        -------
+        d_projectors : list
+            List of <psydac.fem.basic.FemLinearOperator> or <psydac.linalg.basic.LinearOperator>
+            The Dirichlet boundary projectors of each space and in desired form.
+
+        Notes
+        -----
+        See examples/vector_potential_3d.py for a use case of these operators in LinearOperator form.
+
+        """
+        assert kind in ('femlinop', 'linop')
+
+        from psydac.fem.projectors import DirichletMultipatchBoundaryProjector
+        d_projectors = [DirichletMultipatchBoundaryProjector(Vh) for Vh in self.spaces[:-1]]
+
+        if kind == 'femlinop':
+            d_projectors = [FemLinearOperator(fem_domain=Vh, fem_codomain=Vh, linop=d_projector) for Vh, d_projector in zip(self.spaces[:-1], d_projectors)]
+
+        return d_projectors

psydac/feec/conforming_projectors.py

@@ -1211,8 +1211,8 @@ def construct_h1_singlepatch_conforming_projection(Vh, reg_orders=0, p_moments=-
 
 
     def get_vertex_index(coords):
-        nbasis0 = Vh.spaces[coords[0]].nbasis - 1
-        nbasis1 = Vh.spaces[coords[1]].nbasis - 1
+        nbasis0 = Vh.spaces[0].nbasis - 1
+        nbasis1 = Vh.spaces[1].nbasis - 1
 
         # patch local index
         multi_index = [None] * ndim
@@ -1233,6 +1233,9 @@ def vertex_moment_indices(axis, coords, p_moments):
 
     for  co in [(0,0), (1,0), (0,1), (1,1)]:
 
+        if all(Vh.periodic):
+            break
+
         # global index
         ig = get_vertex_index(co)
 
@@ -1304,6 +1307,10 @@ def get_mu_minus(j, coarse_space, fine_space, R):
 
     # boundary condition
     for bn in domain.boundary:
+
+        if Vh.periodic[bn.axis]:
+            continue
+
         space_k = Vh
         axis = bn.axis
 
@@ -1402,6 +1409,9 @@ def edge_moment_index(p, i, axis, ext):
     # boundary condition
     for bn in domain.boundary:
 
+        if Vh.periodic[bn.axis]:
+            continue
+
         axis = bn.axis
         d = 1 - axis
         ext = bn.ext