diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 00a11ea1d..a9746bc34 100644 --- a/psydac/api/ast/utilities.py +++ b/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 diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index e097d4621..0d74f3178 100644 --- a/psydac/api/discretization.py +++ b/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'): """ diff --git a/psydac/api/feec.py b/psydac/api/feec.py index 94f224109..07acc0b4d 100644 --- a/psydac/api/feec.py +++ b/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,6 +156,9 @@ 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)) @@ -137,18 +166,37 @@ def projectors(self, *, kind='global', nquads=None): 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 + diff --git a/psydac/api/fem.py b/psydac/api/fem.py index 098b80904..c96bc9b56 100644 --- a/psydac/api/fem.py +++ b/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: diff --git a/psydac/feec/global_projectors.py b/psydac/feec/global_projectors.py index 773df5cf3..7ecb67a7d 100644 --- a/psydac/feec/global_projectors.py +++ b/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]) + diff --git a/psydac/feec/pull_push.py b/psydac/feec/pull_push.py index c73d9394c..dba0391f5 100644 --- a/psydac/feec/pull_push.py +++ b/psydac/feec/pull_push.py @@ -8,11 +8,12 @@ # ------------------- 'pull_1d_h1', 'pull_1d_l2', + 'pull_2d_h1vec', 'pull_2d_h1', 'pull_2d_hcurl', 'pull_2d_hdiv', 'pull_2d_l2', - 'pull_3d_v', # NOTE: what is this used for? + 'pull_3d_h1vec', # NOTE: what is this used for? 'pull_3d_h1', 'pull_3d_hcurl', 'pull_3d_hdiv', @@ -64,6 +65,35 @@ def f_logical(eta1): #============================================================================== # 2D PULL-BACKS #============================================================================== +def pull_2d_h1vec(f, F): + + assert isinstance(F, BasicCallableMapping) + assert F.ldim == 2 + + f1, f2 = f + + def f1_logical(eta1, eta2): + x, y = F(eta1, eta2) + + a1_phys = f1(x, y) + 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 + return value_1 + + def f2_logical(eta1, eta2): + x, y = F(eta1, eta2) + + a1_phys = f1(x, y) + 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 + return value_2 + + return f1_logical, f2_logical + def pull_2d_h1(f, F): assert isinstance(F, BasicCallableMapping) @@ -161,29 +191,47 @@ def f_logical(eta1, eta2): # TODO [YG 05.10.2022]: # Remove? But it makes sense to return a vector-valued function... -def pull_3d_v(funcs_ini, mapping): +def pull_3d_h1vec(f, F): + + assert isinstance(F, BasicCallableMapping) + assert F.ldim == 3 - mapping = mapping.get_callable_mapping() - f1,f2,f3 = mapping._func_eval - J_inv = mapping._jacobian_inv + f1, f2, f3 = f - def fun(xi1, xi2, xi3): - x = f1(xi1, xi2, xi3) - y = f2(xi1, xi2, xi3) - z = f3(xi1, xi2, xi3) + def f1_logical(eta1, eta2, eta3): + x, y, z = F(eta1, eta2, eta3) - a1_phys = funcs_ini[0](x, y, z) - a2_phys = funcs_ini[1](x, y, z) - a3_phys = funcs_ini[2](x, y, z) + a1_phys = f1(x, y, z) + a2_phys = f2(x, y, z) + a3_phys = f3(x, y, z) - J_inv_value = J_inv(xi1, xi2, xi3) + 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 + return value_1 + + def f2_logical(eta1, eta2, eta3): + x, y, z = F(eta1, eta2, eta3) + + a1_phys = f1(x, y, z) + a2_phys = f2(x, y, z) + 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_3 = J_inv_value[2, 0] * a1_phys + J_inv_value[2, 1] * a2_phys + J_inv_value[2, 2] * a3_phys + return value_2 + + def f3_logical(eta1, eta2, eta3): + x, y, z = F(eta1, eta2, eta3) - return value_1, value_2, value_3 + a1_phys = f1(x, y, z) + a2_phys = f2(x, y, z) + a3_phys = f3(x, y, z) - return fun + 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 + return value_2 + + return f1_logical, f2_logical, f3_logical #============================================================================== def pull_3d_h1(f, F): diff --git a/psydac/feec/tests/test_global_projectors.py b/psydac/feec/tests/test_global_projectors.py index f3b96bd88..f7a3c0cf0 100644 --- a/psydac/feec/tests/test_global_projectors.py +++ b/psydac/feec/tests/test_global_projectors.py @@ -7,6 +7,10 @@ from psydac.fem.tensor import TensorFemSpace from psydac.feec.global_projectors import Projector_H1, Projector_L2 from psydac.ddm.cart import DomainDecomposition +from sympde.topology import Square, Cube +from psydac.api.discretization import discretize +from sympde.topology import element_of, Derham + #============================================================================== @pytest.mark.parametrize('domain', [(0, 2*np.pi)]) @@ -42,7 +46,7 @@ def test_H1_projector_1d(domain, ncells, degree, periodic): # Test if max-norm of error is <= TOL maxnorm_error = abs(vals_u0 - vals_f).max() print(ncells, maxnorm_error) -# assert maxnorm_error <= 1e-14 + assert maxnorm_error <= 1e-9 #============================================================================== @pytest.mark.parametrize('domain', [(0, 2*np.pi)]) @@ -82,7 +86,206 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads): # Test if max-norm of error is <= TOL maxnorm_error = abs(vals_u1 - vals_f).max() print(ncells, maxnorm_error) -# assert maxnorm_error <= 1e-14 + assert maxnorm_error <= 1e-3 + +#============================================================================== +@pytest.mark.parametrize('ncells', [[200,200]]) +@pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]]) +@pytest.mark.parametrize('periodic', [[False, False], [True, True]]) + +def test_derham_projector_2d_hdiv(ncells, degree, periodic): + + domain = Square('Omega', bounds1 = (0,2*np.pi), bounds2 = (0,2*np.pi)) + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + + derham = Derham(domain, ["H1", "Hdiv", "L2"]) + derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True) + P0, P1, P2, PX = derham_h.projectors() + + # Projector onto H1 space (1D interpolation) + + # Function to project + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) + f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) + + # Compute the projection + u0 = P0(f1) + u2 = P2(f1) + u1 = P1((f1,f2)) + ux = PX((f1,f2)) + + # Create evaluation grid, and check if u0(x) == f(x) + xgrid = np.linspace(0, 2*np.pi, num=51) + vals_u0 = np.array([[u0(x, y) for x in xgrid] for y in xgrid]) + vals_u1_1 = np.array([[u1(x, y)[0] for x in xgrid] for y in xgrid]) + vals_u2 = np.array([[u2(x, y) for x in xgrid] for y in xgrid]) + vals_ux_1 = np.array([[ux(x, y)[0] for x in xgrid] for y in xgrid]) + vals_f = np.array([[f1(x, y) for x in xgrid] for y in xgrid]) + + # Test if max-norm of error is <= TOL + maxnorm_error = abs(vals_u0 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_u1_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_u2 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_ux_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + +#============================================================================== +@pytest.mark.parametrize('ncells', [[200,200]]) +@pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]]) +@pytest.mark.parametrize('periodic', [[False, False], [True, True]]) + +def test_derham_projector_2d_hdiv_2(ncells, degree, periodic): + + domain = Square('Omega', bounds1 = (0,1), bounds2 = (0,1)) + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + + derham = Derham(domain, ["H1", "Hdiv", "L2"]) + derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True) + P0, P1, P2, PX = derham_h.projectors() + + # Projector onto H1 space (1D interpolation) + + # Function to project + f1 = lambda xi1, xi2 : xi1**2*(xi1-1.)**2 + #function C0 restricted to [0,1] with periodic BC (0 at x1=0 and x1=1) + f2 = lambda xi1, xi2 : xi2**2*(xi2-1.)**2 + + # Compute the projection + u0 = P0(f1) + u2 = P2(f1) + u1 = P1((f1,f2)) + ux = PX((f1,f2)) + + # Create evaluation grid, and check if u0(x) == f(x) + xgrid = np.linspace(0, 1, num=51) + vals_u0 = np.array([[u0(x, y) for x in xgrid] for y in xgrid]) + vals_u1_1 = np.array([[u1(x, y)[0] for x in xgrid] for y in xgrid]) + vals_u2 = np.array([[u2(x, y) for x in xgrid] for y in xgrid]) + vals_ux_1 = np.array([[ux(x, y)[0] for x in xgrid] for y in xgrid]) + vals_f = np.array([[f1(x, y) for x in xgrid] for y in xgrid]) + + # Test if max-norm of error is <= TOL + maxnorm_error = abs(vals_u0 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_u1_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_u2 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_ux_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + +#============================================================================== +@pytest.mark.parametrize('ncells', [[200,200]]) +@pytest.mark.parametrize('degree', [[2,2], [2,3], [3,3]]) +@pytest.mark.parametrize('periodic', [[False, False], [True, False] ,[True, True]]) + +def test_derham_projector_2d_hcurl(ncells, degree, periodic): + + domain = Square('Omega', bounds1 = (0,2*np.pi), bounds2 = (0,2*np.pi)) + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + + derham = Derham(domain, ["H1", "Hcurl", "L2"]) + derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True) + P0, P1, P2, PX = derham_h.projectors() + + # Projector onto H1 space (1D interpolation) + + # Function to project + f1 = lambda xi1, xi2 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) + f2 = lambda xi1, xi2 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) + + # Compute the projection + u0 = P0(f1) + u2 = P2(f1) + u1 = P1((f1,f2)) + ux = PX((f1,f2)) + + # Create evaluation grid, and check if u0(x) == f(x) + xgrid = np.linspace(0, 2*np.pi, num=51) + vals_u0 = np.array([[u0(x, y) for x in xgrid] for y in xgrid]) + vals_u1_1 = np.array([[u1(x, y)[0] for x in xgrid] for y in xgrid]) + vals_u2 = np.array([[u2(x, y) for x in xgrid] for y in xgrid]) + vals_ux_1 = np.array([[ux(x, y)[0] for x in xgrid] for y in xgrid]) + vals_f = np.array([[f1(x, y) for x in xgrid] for y in xgrid]) + + # Test if max-norm of error is <= TOL + maxnorm_error = abs(vals_u0 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_u1_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_u2 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + maxnorm_error = abs(vals_ux_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 1e-3 + +#============================================================================== +@pytest.mark.parametrize('ncells', [[30,30,30]]) +@pytest.mark.parametrize('degree', [[2,2,2], [2,3,2], [3,3,3]]) +@pytest.mark.parametrize('periodic', [[False, False, False], [True, True, True]]) + +def test_derham_projector_3d(ncells, degree, periodic): + + domain = Cube('Omega', bounds1 = (0,2*np.pi), bounds2 = (0,2*np.pi), bounds3 = (0,2*np.pi)) + domain_h = discretize(domain, ncells=ncells, periodic=periodic) + + derham = Derham(domain) + derham_h = discretize(derham, domain_h, degree=degree, get_H1vec_space = True) + P0, P1, P2, P3, PX = derham_h.projectors() + + # Projector onto H1 space (1D interpolation) + + # Function to project + f1 = lambda xi1, xi2, xi3 : np.sin( xi1 + 0.5 ) * np.cos( xi2 + 0.3 ) * np.sin( 2 * xi3 ) + f2 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.5 ) * np.sin( xi2 - 0.2 ) * np.cos( xi3 ) + f3 = lambda xi1, xi2, xi3 : np.cos( xi1 + 0.7 ) * np.sin( 2*xi2 - 0.2 ) * np.cos( xi3 ) + + # Compute the projection + u0 = P0(f1) + u3 = P3(f1) + u1 = P1((f1,f2,f3)) + u2 = P2((f1,f2,f3)) + ux = PX((f1,f2,f3)) + + # Create evaluation grid, and check if u0(x) == f(x) + xgrid = np.linspace(0, 2*np.pi, num=21) + vals_u0 = np.array([[[u0(x, y, z) for x in xgrid] for y in xgrid] for z in xgrid]) + vals_u1_1 = np.array([[[u1(x, y, z)[0] for x in xgrid] for y in xgrid] for z in xgrid]) + vals_u2_1 = np.array([[[u2(x, y, z)[0] for x in xgrid] for y in xgrid] for z in xgrid]) + vals_ux_1 = np.array([[[ux(x, y, z)[0] for x in xgrid] for y in xgrid] for z in xgrid]) + vals_u3 = np.array([[[u3(x, y, z) for x in xgrid] for y in xgrid] for z in xgrid]) + vals_f = np.array([[[f1(x, y, z) for x in xgrid] for y in xgrid] for z in xgrid]) + + # Test if max-norm of error is <= TOL + maxnorm_error = abs(vals_u0 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 2e-2 + maxnorm_error = abs(vals_u1_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 2e-2 + maxnorm_error = abs(vals_u2_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 2e-2 + maxnorm_error = abs(vals_u3 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 2e-2 + maxnorm_error = abs(vals_ux_1 - vals_f).max() + print(ncells, maxnorm_error) + assert maxnorm_error <= 2e-2 #============================================================================== if __name__ == '__main__': @@ -98,3 +301,14 @@ def test_L2_projector_1d(domain, ncells, degree, periodic, nquads): nquads = degree for nc in ncells: test_L2_projector_1d(domain, nc, degree, periodic, nquads) + + for nc in ncells: + test_derham_projector_2d_hdiv_2([nc, nc], [degree, degree], [periodic, periodic]) + test_derham_projector_2d_hdiv([nc, nc], [degree, degree], [periodic, periodic]) + + for nc in ncells : + test_derham_projector_2d_hcurl([nc, nc], [degree, degree], [periodic, periodic]) + + for nc in ncells[:3] : + test_derham_projector_3d([nc, nc, nc], [degree, degree, degree], [periodic, periodic, periodic]) + diff --git a/psydac/fem/vector.py b/psydac/fem/vector.py index 4ecc896a6..bca983331 100644 --- a/psydac/fem/vector.py +++ b/psydac/fem/vector.py @@ -58,7 +58,8 @@ def __init__( self, *spaces ): self._symbolic_space = None if all(s.symbolic_space for s in spaces): - self._symbolic_space = reduce(lambda x,y:x.symbolic_space*y.symbolic_space, spaces) + symbolic_spaces = [s.symbolic_space for s in spaces] + self._symbolic_space = reduce(lambda x,y:x*y, symbolic_spaces) self._vector_space = BlockVectorSpace(*[V.vector_space for V in self.spaces]) self._refined_space = {} diff --git a/psydac/utilities/utils.py b/psydac/utilities/utils.py index 22179de76..623120569 100644 --- a/psydac/utilities/utils.py +++ b/psydac/utilities/utils.py @@ -59,6 +59,17 @@ def unroll_edges(domain, xgrid): elif xgrid[-1] != xB: return np.array([*xgrid, xgrid[0] + (xB-xA)]) + +#=============================================================================== +def roll_edges(domain, points): + """If necessary, "roll" back intervals that cross boundary of periodic domain. + Changes are made in place to avoid duplicating the array + """ + xA, xB = domain + assert xA < xB + points -=xA + points %=(xB-xA) + points +=xA #=============================================================================== def split_space(Xh): """Split the flattened fem spaces into