Add Hvec space to de Rham sequence, if requested (#282)

The "Hvec" space is the (H^1)^d space taken as the product of d
identical Vh^0 spaces.

List of modifications
---------------------
- Add the "Hvec" space to `DiscreteDerham` objects, if requested to the
`discretize_derham` function through the optional argument
`get_H1vec_space` (default is `False`).
- Add a geometric projection for this space, which is a component-wise
interpolation combined with the pull_back for vector fields.
- Add documentation to the `discretize_derham` function and the
`DiscreteDerham` class.
- Add tests for the subclasses of `GlobalProjector`.

psydac/api/ast/utilities.py

@@ -331,14 +331,14 @@ def compute_atoms_expr_field(atomic_exprs, indices_quad,
         orders = [*get_index(atom).values()]
         args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
         inits += [Assign(test_fun, Mul(*args))]
-        # ...
+        # ...
 
-        # ...
+        # ...
         args     = [IndexedBase(field_name)[idxs], test_fun]
         val_name = SymbolicExpr(atom).name + '_values'
         val      = IndexedBase(val_name)[indices_quad]
         updates += [AugAssign(val,'+',Mul(*args))]
-        # ...
+        # ...
 
     return inits, updates, map_stmts, new_atoms
 
@@ -386,23 +386,23 @@ def compute_atoms_expr_mapping(atomic_exprs, indices_quad,
         element = get_atom_logical_derivatives(atom)
         element_name = 'coeff_' + SymbolicExpr(element).name
 
-        # ...
+        # ...
         test_fun = atom.subs(element, test_function)
         test_fun = SymbolicExpr(test_fun)
-        # ...
+        # ...
 
-        # ...
+        # ...
         orders = [*get_index_logical_derivatives(atom).values()]
         args   = [b[i, d, q] for b, i, d, q in zip(basis, idxs, orders, indices_quad)]
         inits += [Assign(test_fun, Mul(*args))]
-        # ...
+        # ...
 
-        # ...
+        # ...
         val_name = SymbolicExpr(atom).name + '_values'
         val      = IndexedBase(val_name)[indices_quad]
         expr     = IndexedBase(element_name)[idxs] * test_fun
         updates += [AugAssign(val, '+', expr)]
-        # ...
+        # ...
 
     return inits, updates
 

psydac/api/discretization.py

@@ -81,18 +81,42 @@ def change_dtype(V, dtype):
 
     return V
 
-#==============================================================================
-def discretize_derham(derham, domain_h, *args, **kwargs):
+#==============================================================================           
+def discretize_derham(derham, domain_h, get_H1vec_space = False, *args, **kwargs):
+    """
+    Create a discrete De Rham sequence by creating the spaces and then initiating DiscreteDerham object.
+
+    Parameters
+    ----------
+
+    derham : sympde.topology.space.Derham
+        The symbolic Derham sequence
+
+    domain_h   : Geometry
+        Discrete domain where the spaces will be discretized
+        
+    get_H1vec_space : Bool
+        True to also get the "Hvec" space discretizing (H1)^n vector fields
+        
+    **kwargs : list
+        optional parameters for the space discretization
+    """
 
     ldim    = derham.shape
     mapping = domain_h.domain.mapping # NOTE: assuming single-patch domain!
-
     bases  = ['B'] + ldim * ['M']
     spaces = [discretize_space(V, domain_h, basis=basis, **kwargs) \
             for V, basis in zip(derham.spaces, bases)]
 
-    return DiscreteDerham(mapping, *spaces)
+    if get_H1vec_space:
+        X = VectorFunctionSpace('X', domain_h.domain, kind='h1')
+        V0h = spaces[0]
+        Xh  = VectorFemSpace(*([V0h]*ldim))
+        Xh.symbolic_space = X
+        #We still need to specify the symbolic space because of "_recursive_element_of" not implemented in sympde
+        spaces.append(Xh)
 
+    return DiscreteDerham(mapping, *spaces)
 #==============================================================================
 def reduce_space_degrees(V, Vh, *, basis='B', sequence='DR'):
     """

psydac/api/feec.py

@@ -4,25 +4,46 @@
 from psydac.feec.derivatives       import Derivative_1D, Gradient_2D, Gradient_3D
 from psydac.feec.derivatives       import ScalarCurl_2D, VectorCurl_2D, Curl_3D
 from psydac.feec.derivatives       import Divergence_2D, Divergence_3D
-from psydac.feec.global_projectors import Projector_H1, Projector_Hcurl
+from psydac.feec.global_projectors import Projector_H1, Projector_Hcurl, Projector_H1vec
 from psydac.feec.global_projectors import Projector_Hdiv, Projector_L2
 from psydac.feec.pull_push         import pull_1d_h1, pull_1d_l2
-from psydac.feec.pull_push         import pull_2d_h1, pull_2d_hcurl, pull_2d_hdiv, pull_2d_l2
-from psydac.feec.pull_push         import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2
+from psydac.feec.pull_push         import pull_2d_h1, pull_2d_hcurl, pull_2d_hdiv, pull_2d_l2, pull_2d_h1vec
+from psydac.feec.pull_push         import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2, pull_3d_h1vec
+from psydac.fem.vector             import VectorFemSpace
+
 
 __all__ = ('DiscreteDerham',)
 
 #==============================================================================
 class DiscreteDerham(BasicDiscrete):
     """ Represent the discrete De Rham sequence.
+    Should be initialized via discretize_derham function in api.discretization.py
+    
+    Parameters
+    ----------
+
+    mapping : Mapping
+        The mapping from the logical space to the physical space of the discrete De Rham.
+        
+    *spaces : list of FemSpace
+        The discrete spaces of the De Rham sequence
     """
     def __init__(self, mapping, *spaces):
 
         assert (mapping is None) or isinstance(mapping, Mapping)
+
+        self.has_vec = isinstance(spaces[-1], VectorFemSpace)
+
+        if self.has_vec : 
+            dim          = len(spaces) - 2
+            self._spaces = spaces[:-1]
+            self._H1vec    = spaces[-1]
+
+        else :
+            dim           = len(spaces) - 1
+            self._spaces  = spaces
 
-        dim           = len(spaces) - 1
         self._dim     = dim
-        self._spaces  = spaces
         self._mapping = mapping
         self._callable_mapping = mapping.get_callable_mapping() if mapping else None
 
@@ -83,6 +104,11 @@ def V2(self):
     def V3(self):
         return self._spaces[3]
 
+    @property
+    def H1vec(self):
+        assert self.has_vec
+        return self._H1vec
+
     @property
     def spaces(self):
         return self._spaces
@@ -130,25 +156,47 @@ def projectors(self, *, kind='global', nquads=None):
             else:
                 raise TypeError('projector of space type {} is not available'.format(kind))
 
+            if self.has_vec : 
+                Pvec = Projector_H1vec(self.H1vec, nquads)
+
             if self.mapping:
                 P0_m = lambda f: P0(pull_2d_h1(f, self.callable_mapping))
                 P2_m = lambda f: P2(pull_2d_l2(f, self.callable_mapping))
                 if kind == 'hcurl':
                     P1_m = lambda f: P1(pull_2d_hcurl(f, self.callable_mapping))
                 elif kind == 'hdiv':
                     P1_m = lambda f: P1(pull_2d_hdiv(f, self.callable_mapping))
-                return P0_m, P1_m, P2_m
-            return P0, P1, P2
+                if self.has_vec : 
+                    Pvec_m = lambda f: Pvec(pull_2d_h1vec(f, self.callable_mapping))
+                    return P0_m, P1_m, P2_m, Pvec_m
+                else : 
+                    return P0_m, P1_m, P2_m
+
+            if self.has_vec :
+                return P0, P1, P2, Pvec
+            else : 
+                return P0, P1, P2
 
         elif self.dim == 3:
             P0 = Projector_H1   (self.V0)
             P1 = Projector_Hcurl(self.V1, nquads)
             P2 = Projector_Hdiv (self.V2, nquads)
             P3 = Projector_L2   (self.V3, nquads)
+            if self.has_vec : 
+                Pvec = Projector_H1vec(self.H1vec)
             if self.mapping:
                 P0_m = lambda f: P0(pull_3d_h1   (f, self.callable_mapping))
                 P1_m = lambda f: P1(pull_3d_hcurl(f, self.callable_mapping))
                 P2_m = lambda f: P2(pull_3d_hdiv (f, self.callable_mapping))
                 P3_m = lambda f: P3(pull_3d_l2   (f, self.callable_mapping))
-                return P0_m, P1_m, P2_m, P3_m
-            return P0, P1, P2, P3
+                if self.has_vec : 
+                    Pvec_m = lambda f: Pvec(pull_3d_h1vec(f, self.callable_mapping))
+                    return P0_m, P1_m, P2_m, P3_m, Pvec_m
+                else : 
+                    return P0_m, P1_m, P2_m, P3_m
+
+            if self.has_vec :
+                return P0, P1, P2, P3, Pvec
+            else : 
+                return P0, P1, P2, P3
+

psydac/api/fem.py

@@ -1539,7 +1539,7 @@ def assemble(self, **kwargs):
             v = kwargs[key]
             if isinstance(v, FemField):
                 if not v.coeffs.ghost_regions_in_sync:
-                        v.coeffs.update_ghost_regions()
+                    v.coeffs.update_ghost_regions()
                 if v.space.is_product:
                     coeffs = v.coeffs
                     if self._symbolic_space.is_broken:

psydac/feec/global_projectors.py

@@ -12,6 +12,8 @@
 from psydac.fem.tensor import TensorFemSpace
 from psydac.fem.vector import VectorFemSpace
 
+from psydac.utilities.utils import roll_edges
+
 from abc import ABCMeta, abstractmethod
 
 __all__ = ('GlobalProjector', 'Projector_H1', 'Projector_Hcurl', 'Projector_Hdiv', 'Projector_L2',
@@ -147,6 +149,11 @@ def __init__(self, space, nquads = None):
                     if quad_x[j] is None:
                         u, w = uw[j]
                         global_quad_x, global_quad_w = quadrature_grid(V.histopolation_grid, u, w)
+                        #"roll" back points to the interval to ensure that the quadrature points are
+                        #in the domain. Only usefull in th eperiodic case (else do nothing)
+                        #if not used then you will have quadrature points outside of the domain which 
+                        #might cause problem when your function is only defined inside the domain
+                        roll_edges(V.domain, global_quad_x) 
                         quad_x[j] = global_quad_x[s:e+1]
                         quad_w[j] = global_quad_w[s:e+1]
                     local_x, local_w = quad_x[j], quad_w[j]
@@ -564,6 +571,70 @@ def __call__(self, fun):
         """
         return super().__call__(fun)
 
+class Projector_H1vec(GlobalProjector):
+    """
+    Projector from H1^3 = H1 x H1 x H1 to a conforming finite element space, i.e.
+    a finite dimensional subspace of H1^3, constructed with tensor-product
+    B-splines in 2 or 3 dimensions.
+    This is a global projector constructed over a tensor-product grid in the
+    logical domain. The vertices of this grid are obtained as the tensor
+    product of the 1D splines' Greville points along each direction.
+    
+    Parameters
+    ----------
+    H1vec : ProductFemSpace
+        H1 x H1 x H1-conforming finite element space, codomain of the projection
+        operator.
+        
+    nquads : list(int) | tuple(int)
+        Number of quadrature points along each direction, to be used in Gauss
+        quadrature rule for computing the (approximated) degrees of freedom.
+    """
+    def _structure(self, dim):
+        if dim == 3:
+            return [
+                ['I', 'I', 'I'],
+                ['I', 'I', 'I'],
+                ['I', 'I', 'I']
+            ]
+        elif dim == 2:
+            return [
+                ['I', 'I'],
+                ['I', 'I']
+            ]
+        else:
+            raise NotImplementedError('The H1vec projector is only available in 2D or 3D.')
+
+    def _function(self, dim):
+        if dim == 3: return evaluate_dofs_3d_vec
+        elif dim == 2: return evaluate_dofs_2d_vec
+        else:
+            raise NotImplementedError('The H1vec projector is only available in 2/3D.')
+
+    #--------------------------------------------------------------------------
+    def __call__(self, fun):
+        r"""
+        Project vector function onto the H1 x H1 x H1-conforming finite element
+        space. This happens in the logical domain $\hat{\Omega}$.
+        
+        Parameters
+        ----------
+        fun : list/tuple of callables
+            Scalar components of the real-valued vector function to be
+            projected, with arguments the coordinates (x_1, ..., x_N) of a
+            point in the logical domain. These correspond to the coefficients
+            of a vector-field.
+            $fun_i : \hat{\Omega} \mapsto \mathbb{R}$ with i = 1, ..., N.
+            
+        Returns
+        -------
+        field : FemField
+            Field obtained by projection (element of the H1^3-conforming
+            finite element space). This is also a real-valued vector function
+            in the logical domain.
+        """
+        return super().__call__(fun)
+
 #==============================================================================
 # 1D DEGREES OF FREEDOM
 #==============================================================================
@@ -673,6 +744,24 @@ def evaluate_dofs_2d_2form(
                     F[i1, i2] += quad_w1[i1, g1] * quad_w2[i2, g2] * \
                             f(quad_x1[i1, g1], quad_x2[i2, g2])
 
+#------------------------------------------------------------------------------    
+def evaluate_dofs_2d_vec(
+        intp_x1, intp_x2,      # interpolation points
+        F1, F2,                # array of degrees of freedom (intent out)
+        f1, f2,                # input scalar function (callable)
+        ):
+
+    n1, n2 = F1.shape
+    for i1 in range(n1):
+        for i2 in range(n2):
+            F1[i1, i2] = f1(intp_x1[i1], intp_x2[i2])
+
+    n1, n2 = F2.shape
+    for i1 in range(n1):
+        for i2 in range(n2):
+            F2[i1, i2] = f2(intp_x1[i1], intp_x2[i2])
+
+
 #==============================================================================
 # 3D DEGREES OF FREEDOM
 #==============================================================================
@@ -791,3 +880,30 @@ def evaluate_dofs_3d_3form(
                             F[i1, i2, i3] += \
                                     quad_w1[i1, g1] * quad_w2[i2, g2] * quad_w3[i3, g3] * \
                                     f(quad_x1[i1, g1], quad_x2[i2, g2], quad_x3[i3, g3])
+
+#------------------------------------------------------------------------------
+def evaluate_dofs_3d_vec(
+        intp_x1, intp_x2, intp_x3, # interpolation points
+        F1, F2, F3,                # array of degrees of freedom (intent out)
+        f1, f2, f3,                # input scalar function (callable)
+        ):
+
+    # evaluate input functions at interpolation points (make sure that points are in [0, 1])
+    n1, n2, n3 = F1.shape
+    for i1 in range(n1):
+        for i2 in range(n2):
+            for i3 in range(n3):
+                F1[i1, i2, i3] = f1(intp_x1[i1], intp_x2[i2], intp_x3[i3])
+
+    n1, n2, n3 = F2.shape
+    for i1 in range(n1):
+        for i2 in range(n2):
+            for i3 in range(n3):
+                F2[i1, i2, i3] = f2(intp_x1[i1], intp_x2[i2], intp_x3[i3])
+
+    n1, n2, n3 = F3.shape
+    for i1 in range(n1):
+        for i2 in range(n2):
+            for i3 in range(n3):
+                F3[i1, i2, i3] = f3(intp_x1[i1], intp_x2[i2], intp_x3[i3])
+