answering Yaman review

psydac/api/discretization.py

@@ -116,7 +116,7 @@ def discretize_derham(derham, domain_h, get_H1vec_space = False, *args, **kwargs
         #We still need to specify the symbolic space because of "_recursive_element_of" not implemented in sympde
         spaces.append(Xh)
 
-    return DiscreteDerham(mapping, get_H1vec_space, *spaces)
+    return DiscreteDerham(mapping, *spaces)
 #==============================================================================
 def reduce_space_degrees(V, Vh, *, basis='B', sequence='DR'):
     """

psydac/api/feec.py

@@ -9,6 +9,8 @@
 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, pull_2d_vec
 from psydac.feec.pull_push         import pull_3d_h1, pull_3d_hcurl, pull_3d_hdiv, pull_3d_l2, pull_3d_vec
+from psydac.fem.vector             import VectorFemSpace
+
 
 __all__ = ('DiscreteDerham',)
 
@@ -22,23 +24,20 @@ class DiscreteDerham(BasicDiscrete):
 
     mapping : Mapping
         The mapping from the logical space to the physical space of the discrete De Rham.
-
-    get_vec : Bool
-        True to also get the "Hvec" space discretizing (H1)^n vector fields
         
     *spaces : list of 
         The discrete spaces of the De Rham sequence
     """
-    def __init__(self, mapping, get_H1vec_space=False, *spaces):
+    def __init__(self, mapping, *spaces):
 
         assert (mapping is None) or isinstance(mapping, Mapping)
 
-        self.has_vec = get_H1vec_space
+        self.has_vec = isinstance(spaces[-1], VectorFemSpace)
 
         if self.has_vec : 
             dim          = len(spaces) - 2
             self._spaces = spaces[:-1]
-            self._Vvec    = spaces[-1]
+            self._H1vec    = spaces[-1]
 
         else :
             dim           = len(spaces) - 1
@@ -106,9 +105,9 @@ def V3(self):
         return self._spaces[3]
 
     @property
-    def Vvec(self):
+    def H1vec(self):
         assert self.has_vec
-        return self._Vvec
+        return self._H1vec
 
     @property
     def spaces(self):
@@ -158,7 +157,7 @@ def projectors(self, *, kind='global', nquads=None):
                 raise TypeError('projector of space type {} is not available'.format(kind))
 
             if self.has_vec : 
-                Pvec = Projector_H1vec(self.Vvec, nquads)
+                Pvec = Projector_H1vec(self.H1vec, nquads)
 
             if self.mapping:
                 P0_m = lambda f: P0(pull_2d_h1(f, self.callable_mapping))
@@ -184,7 +183,7 @@ def projectors(self, *, kind='global', nquads=None):
             P2 = Projector_Hdiv (self.V2, nquads)
             P3 = Projector_L2   (self.V3, nquads)
             if self.has_vec : 
-                Pvec = Projector_H1vec(self.Vvec)
+                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))

psydac/feec/global_projectors.py

@@ -577,11 +577,16 @@ class Projector_H1vec(GlobalProjector):
     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:
@@ -609,6 +614,7 @@ 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
@@ -617,6 +623,7 @@ def __call__(self, fun):
             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

psydac/feec/pull_push.py

@@ -79,7 +79,7 @@ def f1_logical(eta1, eta2):
         a2_phys = f2(x, y)
 
         J_inv_value = F.jacobian_inv(eta1, eta2)
-        value_1 = J_inv_value[0,0]*a1_phys + J_inv_value[0,1]*a2_phys
+        value_1 = J_inv_value[0, 0] * a1_phys + J_inv_value[0, 1] * a2_phys
         return value_1
 
     def f2_logical(eta1, eta2):
@@ -89,7 +89,7 @@ def f2_logical(eta1, eta2):
         a2_phys = f2(x, y)
 
         J_inv_value = F.jacobian_inv(eta1, eta2)
-        value_2 = J_inv_value[1,0]*a1_phys + J_inv_value[1,1]*a2_phys
+        value_2 = J_inv_value[1, 0] * a1_phys + J_inv_value[1, 1] * a2_phys
         return value_2
 
     return f1_logical, f2_logical
@@ -206,7 +206,7 @@ def f1_logical(eta1, eta2, eta3):
         a3_phys = f3(x, y, z)
 
         J_inv_value = F.jacobian_inv(eta1, eta2, eta3)
-        value_1 = J_inv_value[0,0]*a1_phys + J_inv_value[0,1]*a2_phys + J_inv_value[0,2]*a3_phys
+        value_1 = J_inv_value[0, 0] * a1_phys + J_inv_value[0, 1] * a2_phys + J_inv_value[0, 2] * a3_phys
         return value_1
 
     def f2_logical(eta1, eta2, eta3):
@@ -217,7 +217,7 @@ def f2_logical(eta1, eta2, eta3):
         a3_phys = f3(x, y, z)
 
         J_inv_value = F.jacobian_inv(eta1, eta2, eta3)
-        value_2 = J_inv_value[1,0]*a1_phys + J_inv_value[1,1]*a2_phys + J_inv_value[1,2]*a3_phys
+        value_2 = J_inv_value[1, 0] * a1_phys + J_inv_value[1, 1] * a2_phys + J_inv_value[1, 2] * a3_phys
         return value_2
 
     def f3_logical(eta1, eta2, eta3):
@@ -228,7 +228,7 @@ def f3_logical(eta1, eta2, eta3):
         a3_phys = f3(x, y, z)
 
         J_inv_value = F.jacobian_inv(eta1, eta2, eta3)
-        value_2 = J_inv_value[2,0]*a1_phys + J_inv_value[2,1]*a2_phys + J_inv_value[2,2]*a3_phys
+        value_2 = J_inv_value[2, 0] * a1_phys + J_inv_value[2, 1] * a2_phys + J_inv_value[2, 2] * a3_phys
         return value_2
 
     return f1_logical, f2_logical, f3_logical