[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Tue, 27 Aug 2019 14:37:04 -0400 (EDT) |
branch: master
commit ea2609cca7cc7387a020685025a2e73be3cae0bb
Author: Yves Renard <address@hidden>
Date: Thu Aug 22 20:41:13 2019 +0200
HHO stabilization operator on a target space different of the original one
---
doc/sphinx/source/userdoc/hho.rst | 15 +-
interface/tests/python/demo_elasticity_HHO.py | 48 +++--
interface/tests/python/demo_laplacian_HHO.py | 43 ++--
src/getfem_HHO.cc | 275 +++++++++++++++-----------
src/getfem_fem.cc | 158 +++++++++++++--
5 files changed, 381 insertions(+), 158 deletions(-)
diff --git a/doc/sphinx/source/userdoc/hho.rst
b/doc/sphinx/source/userdoc/hho.rst
index a5a3086..42e47e9 100644
--- a/doc/sphinx/source/userdoc/hho.rst
+++ b/doc/sphinx/source/userdoc/hho.rst
@@ -20,11 +20,16 @@ HHO elements
HHO elements are composite ones having a polynomial approximation space for
the interior of the element and a polynomial approximation for each face of the
element. Moreover, this is a discontinous approximation, in the sens that no
continuity is prescribed between the approximation inside the element and the
approximation on the faces, neither than between the approximations on two
different faces of the element. However, when two neighbour elements share a
face, the approximation on th [...]
- getfem::pfem pf = getfem::fem_descriptor("HHO(FEM_DISCONTINOUS_PK(2,2,0.1),
FEM_DISCONTINOUS_PK(1,2,0.1))");
+ getfem::pfem pf = getfem::fem_descriptor("HHO(FEM_SIMPLEX_IPK(2,2),
FEM_SIMPLEX_CIPK(1,2))");
-The first argument to ``FEM_HHO(...)`` is the fem for the interior of the
element. It has to be a discontinuous FEM. The method
``FEM_DISCONTINOUS_PK(2,2,0.1)`` is a discontinous method having its degrees of
freedom in the strict interior of the element, which ensure that no dof
identification will be done. The second argument is the fem for the faces (if
only one method is given, it will be applied to all faces, but it is also
possible to give a different method for each face). Their is [...]
+The first argument to ``FEM_HHO(...)`` is the fem for the interior of the
element. It has to be a discontinuous FEM. The method ``FEM_SIMPLEX_IPK(2,2)``
is a discontinous method having its degrees of freedom in the strict interior
of the element, which ensure that no dof identification will be done. The
second argument is the fem for the faces (if only one method is given, it will
be applied to all faces, but it is also possible to give a different method for
each face). Their is no veri [...]
-+ ``FEM_IPK`` and PK on quadrilateral ...
+For the moment, the fursnished element for interior and faces are
+- ``FEM_SIMPLEX_IPK(n,k)`` : interior PK element of degree k for the simplices
in dimension n (equivalent to ``FEM_PK_DISCONTINUOUS(n,k,0.1)``).
+- ``FEM_QUAD_IPK(n,k)`` : interior PK element of degree k for the
quadrilaterals in dimension n.
+- ``FEM_PRISM_IPK(n,k)`` : interior PK element of degree k for the prisms in
dimension n.
+- ``FEM_SIMPLEX_CIPK(n,k)`` : interior PK element on simplices which is
additionnaly connectable. Designed to be use on HHO element face.
+- ``FEM_QUAD_CIPK(k)`` : interior PK element on a quadrilateral which is
additionnaly connectable. Designed to be use on HHO element face.
Reconstruction operators
------------------------
@@ -51,7 +56,7 @@ This is an example of use with the Python interface for a
two-dimensional triang
mfu = gf.MeshFem(m, 1)
mfgu = gf.MeshFem(m, N)
-
mfu.set_fem(gf.Fem('FEM_HHO(FEM_PK_DISCONTINUOUS(2,2,0.1),FEM_PK_DISCONTINUOUS(1,2,0.1))'))
+ mfu.set_fem(gf.Fem('FEM_HHO(FEM_SIMPLEX_IPK(2,2),FEM_SIMPLEX_CIPK(1,2))'))
mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
md = gf.Model('real')
@@ -149,4 +154,4 @@ The corresponding elementary transformations can be added
to the model by the tw
add_HHO_stabilization(model, transname);
add_HHO_symmetrized_stabilization(model, transname);
-and used into the weak form language as ``Elementary_transformation(u,
HHO_stab)``, if ``transname="HHO_stab"``. There is no third argument for these
transformations since the arrival space is the one of the variable. An example
of use is also given in the test programs
:file:`interface/tests/demo_laplacian_HHO.py` and
:file:`interface/tests/demo_elasticity_HHO.py`.
+and used into the weak form language as ``Elementary_transformation(u,
HHO_stab)``, if ``transname="HHO_stab"``. A third argument is optional to
specify the target (HHO) space (the default is one of the variable itself). An
example of use is also given in the test programs
:file:`interface/tests/demo_laplacian_HHO.py` and
:file:`interface/tests/demo_elasticity_HHO.py`.
diff --git a/interface/tests/python/demo_elasticity_HHO.py
b/interface/tests/python/demo_elasticity_HHO.py
index c3653d9..fbaae48 100644
--- a/interface/tests/python/demo_elasticity_HHO.py
+++ b/interface/tests/python/demo_elasticity_HHO.py
@@ -35,7 +35,7 @@ NX = 20 # Mesh parameter.
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
-using_HHO = True # Use HHO method or standard Lagrange FEM
+using_HHO = False # Use HHO method or standard Lagrange FEM
using_symmetric_gradient = True # Use symmetric gradient reconstruction or
not
E = 1 # Young's modulus
@@ -43,12 +43,18 @@ nu = 0.499 # Poisson ratio
cmu = E/(2*(1+nu)) # Lame coefficient
clambda = 2*cmu*nu/(1-2*nu) # Lame coefficient
+use_quad = False # Quadrilaterals or triangles
# Create a simple cartesian mesh
-m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
np.arange(0,1+1./NX,1./NX))
-# m = gf.Mesh('cartesian',[0:1/NX:1],[0:1/NX:1]);
-#
m=gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[%d,%d];'
% (NX, NX))
-#
m=gf.Mesh('import','structured','GT="GT_QK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[1,1];')
+if (use_quad):
+ m = gf.Mesh('cartesian', np.arange(0,1+1./NX,1./NX),
+ np.arange(0,1+1./NX,1./NX));
+ #
m=gf.Mesh('import','structured','GT="GT_QK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[1,1];')
+else:
+ m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
+ np.arange(0,1+1./NX,1./NX))
+ #
m=gf.Mesh('import','structured','GT="GT_PK(2,1)";SIZES=[1,1];NOISED=1;NSUBDIV=[%d,%d];'
% (NX, NX))
+
N = m.dim();
@@ -59,19 +65,35 @@ mfur = gf.MeshFem(m, N)
mfrhs = gf.MeshFem(m, 1)
if (using_HHO):
-
mfu.set_fem(gf.Fem('FEM_HHO(FEM_PK_DISCONTINUOUS(2,2,0.1),FEM_PK_DISCONTINUOUS(1,2,0.1))'))
- mfur.set_fem(gf.Fem('FEM_PK(2,3)'))
+ if (use_quad):
+ mfu.set_fem(gf.Fem('FEM_HHO(FEM_QUAD_IPK(2,2),FEM_SIMPLEX_CIPK(1,2))'))
+ #
mfu.set_fem(gf.Fem('FEM_HHO(FEM_QK_DISCONTINUOUS(2,2,0.1),FEM_SIMPLEX_CIPK(1,2))'))
+ mfur.set_fem(gf.Fem('FEM_QK(2,3)'))
+ else:
+ mfu.set_fem(gf.Fem('FEM_HHO(FEM_SIMPLEX_IPK(2,2),FEM_SIMPLEX_CIPK(1,2))'))
+ mfur.set_fem(gf.Fem('FEM_PK(2,3)'))
else:
- mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
- mfur.set_fem(gf.Fem('FEM_PK(2,2)'))
-
-mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
-mfrhs.set_fem(gf.Fem('FEM_PK(2,2)'))
+ if (use_quad):
+ mfu.set_fem(gf.Fem('FEM_QK(2,2)'))
+ mfur.set_fem(gf.Fem('FEM_QK(2,2)'))
+ else:
+ mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
+ mfur.set_fem(gf.Fem('FEM_PK(2,2)'))
+
+if (use_quad):
+ mfgu.set_fem(gf.Fem('FEM_QK(2,2)'))
+ mfrhs.set_fem(gf.Fem('FEM_QK(2,2)'))
+else:
+ mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
+ mfrhs.set_fem(gf.Fem('FEM_PK(2,2)'))
print('nbdof : %d' % mfu.nbdof());
# Integration method used
-mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
+if (use_quad):
+ mim = gf.MeshIm(m, gf.Integ('IM_GAUSS_PARALLELEPIPED(2,4)'))
+else:
+ mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
# Boundary selection
flst = m.outer_faces()
diff --git a/interface/tests/python/demo_laplacian_HHO.py
b/interface/tests/python/demo_laplacian_HHO.py
index db8de14..c6f287a 100644
--- a/interface/tests/python/demo_laplacian_HHO.py
+++ b/interface/tests/python/demo_laplacian_HHO.py
@@ -19,7 +19,7 @@
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
-""" 2D Poisson problem using HHO methods.
+""" Discretization of a Poisson problem using HHO methods.
This program is used to check that python-getfem is working. This is
also a good example of use of GetFEM++.
@@ -31,16 +31,22 @@ import getfem as gf
import numpy as np
## Parameters
-NX = 10 # Mesh parameter.
+NX = 20 # Mesh parameter.
+N = 2
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
using_HHO = True # Use HHO method or standard Lagrange FEM
# Create a simple cartesian mesh
-m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
- np.arange(0,1+1./NX,1./NX))
-N = m.dim();
+if (N == 2):
+ m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
+ np.arange(0,1+1./NX,1./NX))
+elif (N == 3):
+ m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
+ np.arange(0,1+1./NX,1./NX), np.arange(0,1+1./NX,1./NX))
+
+
# Meshfems
mfu = gf.MeshFem(m, 1)
@@ -49,19 +55,22 @@ mfur = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
if (using_HHO):
-
mfu.set_fem(gf.Fem('FEM_HHO(FEM_PK_DISCONTINUOUS(2,2,0.1),FEM_PK_DISCONTINUOUS(1,2,0.1))'))
- mfur.set_fem(gf.Fem('FEM_PK(2,3)'))
+ mfu.set_fem(gf.Fem('FEM_HHO(FEM_SIMPLEX_IPK(%d,2),FEM_SIMPLEX_CIPK(%d,2))' %
(N,N-1)))
+ mfur.set_fem(gf.Fem('FEM_PK(%d,3)' % N))
else:
- mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
- mfur.set_fem(gf.Fem('FEM_PK(2,2)'))
+ mfu.set_fem(gf.Fem('FEM_PK(%d,2)' % N))
+ mfur.set_fem(gf.Fem('FEM_PK(%d,2)' % N))
-mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
-mfrhs.set_fem(gf.Fem('FEM_PK(2,2)'))
+mfgu.set_fem(gf.Fem('FEM_PK(%d,2)' % N))
+mfrhs.set_fem(gf.Fem('FEM_PK(%d,2)' % N))
print('nbdof : %d' % mfu.nbdof());
# Integration method used
-mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
+if (N == 2):
+ mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
+else:
+ mim = gf.MeshIm(m, gf.Integ('IM_TETRAHEDRON(5)'))
# Boundary selection
flst = m.outer_faces()
@@ -90,7 +99,10 @@ Ue = mfur.eval('y*(y-1)*x*(x-1)+x*x*x*x*x')
# Interpolate the source term
F1 = mfrhs.eval('-(2*(x*x+y*y)-2*x-2*y+20*x*x*x)')
-F2 = mfrhs.eval('[y*(y-1)*(2*x-1) + 5*x*x*x*x, x*(x-1)*(2*y-1)]')
+if (N == 2):
+ F2 = mfrhs.eval('[y*(y-1)*(2*x-1) + 5*x*x*x*x, x*(x-1)*(2*y-1)]')
+else:
+ F2 = mfrhs.eval('[y*(y-1)*(2*x-1) + 5*x*x*x*x, x*(x-1)*(2*y-1), 0]')
# Model
md = gf.Model('real')
@@ -131,8 +143,7 @@ md.add_source_term_brick(mim, 'u', 'VolumicData')
# Neumann condition.
md.add_initialized_fem_data('NeumannData', mfrhs, F2)
-md.add_normal_source_term_brick(mim, 'u', 'NeumannData',
- NEUMANN_BOUNDARY_NUM)
+md.add_normal_source_term_brick(mim, 'u', 'NeumannData', NEUMANN_BOUNDARY_NUM)
# Dirichlet condition on the left.
md.add_initialized_fem_data("Ue", mfur, Ue)
@@ -179,6 +190,6 @@ print('You can view the solution with (for example):')
print('gmsh laplacian.pos')
-if (H1error > 1e-3):
+if (H1error > 3e-5):
print('Error too large !')
exit(1)
diff --git a/src/getfem_HHO.cc b/src/getfem_HHO.cc
index 7831047..55266cd 100644
--- a/src/getfem_HHO.cc
+++ b/src/getfem_HHO.cc
@@ -33,6 +33,7 @@ namespace getfem {
// problem can be solved componentwise can be more exploited in
// avoiding the computation of the whole matrix M2.
// - The vectorization can be avoided in most cases
+ // - Avoid some multiple ctx
class _HHO_reconstructed_gradient
@@ -651,7 +652,7 @@ namespace getfem {
size_type cv, base_matrix &M) const {
// The reconstructed variable "S" is described on mf2 and computed by
// S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)))
- // where P__{\dT} et P_T are L2 projections on the boundary and on the
+ // where P_{\dT} et P_T are L2 projections on the boundary and on the
// interior of T on the corresponding discrete spaces.
// Note that P_{\dT}(v_{dT}) = v_{dT} and P_T(v_T) = v_T and D is
// the reconstructed value on P^{k+1} given by the formula:
@@ -664,21 +665,20 @@ namespace getfem {
// "v_{dT}" is the field on the boundary of T and "n" is the outward
// unit normal.
// The implemented formula is
- // S(v) = v_{dT} - P_{\dT}D(v) - P_{\dT}(v_T) + P_{\dT}(P_T(D(v)))
- // by the mean of the projection matrix from P^{k+1} to the original
space
+ // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)) )
+ // by the mean of the projection matrix from P^{k+1} to the target space
// and the projection matrix from interior space to the boundary space.
// As it is built, S(v) is zero on interior dofs.
- GMM_ASSERT1(&mf1 == &mf2, "The HHO stabilization transformation is "
- "only defined on the HHO space to itself");
-
// Obtaining the fem descriptors
pfem pf1 = mf1.fem_of_element(cv);
short_type degree = pf1->estimated_degree();
bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
pfem pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be
// changed for an interior PK method
- pfem pfi = interior_fem_of_hho_method(pf1);
+ pfem pf3 = mf2.fem_of_element(cv);
+ pfem pf1i = interior_fem_of_hho_method(pf1);
+ pfem pf3i = interior_fem_of_hho_method(pf3);
papprox_integration pim
= classical_approx_im(pgt, dim_type(2*degree+2))->approx_method();
@@ -690,11 +690,15 @@ namespace getfem {
= HHO_pgp_pool(pgt, pim->pintegration_points());
pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
- pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
+ pfem_precomp pfp3 = HHO_pfp_pool(pf3, pim->pintegration_points());
+ pfem_precomp pfp1i = HHO_pfp_pool(pf1i, pim->pintegration_points());
+ pfem_precomp pfp3i = HHO_pfp_pool(pf3i, pim->pintegration_points());
fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
- fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
+ fem_interpolation_context ctx3(pgp, pfp3, 0, G, cv);
+ fem_interpolation_context ctx1i(pgp, pfp1i, 0, G, cv);
+ fem_interpolation_context ctx3i(pgp, pfp3i, 0, G, cv);
size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
base_vector un(N);
@@ -702,23 +706,27 @@ namespace getfem {
size_type ndof1 = pf1->nb_dof(cv) * qmult1;
size_type qmult2 = Q / pf2->target_dim();
size_type ndof2 = pf2->nb_dof(cv) * qmult2;
- size_type qmulti = Q / pfi->target_dim();
- size_type ndofi = pfi->nb_dof(cv) * qmulti;
+ size_type qmult3 = Q / pf3->target_dim();
+ size_type ndof3 = pf3->nb_dof(cv) * qmult3;
+ size_type qmult1i = Q / pf1i->target_dim();
+ size_type ndof1i = pf1i->nb_dof(cv) * qmult1i;
+ size_type qmult3i = Q / pf3i->target_dim();
+ size_type ndof3i = pf3i->nb_dof(cv) * qmult3i;
- base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
- base_matrix tv1p, tv2p, tvi;
+ base_tensor t1, t2, t3, t1i, t3i, tv1, tv2, t1p, t2p;
+ base_matrix tv1p, tv2p, tv3p, tv1i, tv3i;
base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
base_matrix M3(Q, ndof1), M4(Q, ndof2);
base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
- base_matrix M7(ndof1, ndof1), M7inv(ndof1, ndof1), M8(ndof1, ndof2);
- base_matrix M9(ndof1, ndof1), MD(ndof2, ndof1);
+ base_matrix M7(ndof3, ndof3), M7inv(ndof3, ndof3), M8(ndof3, ndof2);
+ base_matrix M9(ndof3, ndof1), M10(ndof3, ndof3), MD(ndof2, ndof1);
// Integrals on the element : \int_T Grad(D).Grad(w) (M2)
// \int_T Grad(v_T).Grad(w) (M1)
// \int_T D (M4) and \int_T v_T (M3)
for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
- ctx1.set_ii(ipt); ctx2.set_ii(ipt);
+ ctx1.set_ii(ipt); ctx2.set_ii(ipt); ctx3.set_ii(ipt);
scalar_type coeff = pim->coeff(ipt) * ctx1.J();
ctx1.grad_base_value(t1);
@@ -733,6 +741,9 @@ namespace getfem {
ctx2.base_value(t2p);
vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
+ ctx3.base_value(t3);
+ vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
+
for (size_type i = 0; i < ndof2; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
for (size_type k = 0; k < Q*N; ++k)
@@ -753,26 +764,31 @@ namespace getfem {
for (size_type k = 0; k < Q; ++k)
M3(k, i) += coeff * tv1p(i, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
- for (size_type j = 0; j < ndof1; ++j)
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof3; ++j)
for (size_type k = 0; k < Q; ++k)
- M7(i, j) += coeff * tv1p(i, k) * tv1p(j, k);
+ M7(i, j) += coeff * tv3p(i, k) * tv3p(j, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
for (size_type k = 0; k < Q; ++k)
- M8(i, j) += coeff * tv1p(i, k) * tv2p(j, k);
+ M8(i, j) += coeff * tv3p(i, k) * tv2p(j, k);
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof1; ++j)
+ for (size_type k = 0; k < Q; ++k)
+ M9(i, j) += coeff * tv3p(i, k) * tv1p(j, k);
}
// Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
- ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
+ ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctx3.set_face_num(ifc);
+ ctx1i.set_face_num(ifc); ctx3i.set_face_num(ifc);
size_type first_ind = pim->ind_first_point_on_face(ifc);
for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
- ctx1.set_ii(first_ind+ipt);
- ctx2.set_ii(first_ind+ipt);
- ctxi.set_ii(first_ind+ipt);
+ ctx1.set_ii(first_ind+ipt); ctx2.set_ii(first_ind+ipt);
+ ctx3.set_ii(first_ind+ipt); ctx1i.set_ii(first_ind+ipt);
+ ctx3i.set_ii(first_ind+ipt);
scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
ctx2.grad_base_value(t2);
@@ -784,8 +800,14 @@ namespace getfem {
ctx1.base_value(t1);
vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
- ctxi.base_value(ti);
- vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), Q);
+ ctx1i.base_value(t1i);
+ vectorize_base_tensor(t1i, tv1i, ndof1i, pf1i->target_dim(), Q);
+
+ ctx3i.base_value(t3i);
+ vectorize_base_tensor(t3i, tv3i, ndof3i, pf3i->target_dim(), Q);
+
+ ctx3.base_value(t3);
+ vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
@@ -793,29 +815,34 @@ namespace getfem {
for (size_type j = 0; j < ndof2; ++j)
for (size_type k1 = 0; k1 < Q; ++k1) {
scalar_type b(0), a = coeff *
- (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
+ (tv1p(i, k1) - (i < ndof1i ? tv1i(i, k1) : 0.));
for (size_type k2 = 0; k2 < N; ++k2)
b += a * tv2.as_vector()[j+(k1 + k2*Q)*ndof2] * un[k2];
M1(j, i) += b;
}
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
- for (size_type j = 0; j < ndof1; ++j)
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof3; ++j)
for (size_type k = 0; k < Q; ++k)
- M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
+ M7(i, j) += coeff * tv3p(i,k) * tv3p(j, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
for (size_type k = 0; k < Q; ++k)
- M8(i, j) += coeff * tv1p(i,k) * tv2p(j, k);
+ M8(i, j) += coeff * tv3p(i,k) * tv2p(j, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
- for (size_type j = 0; j < ndofi; ++j)
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof1; ++j)
+ for (size_type k = 0; k < Q; ++k)
+ M9(i, j) += coeff * tv3p(i, k) * tv1p(j, k);
+
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof3i; ++j)
for (size_type k = 0; k < Q; ++k)
- M9(i, j) += coeff * tv1p(i,k) * tvi(j, k);
+ M10(i, j) += coeff * tv3p(i,k) * tv3i(j, k);
}
}
-
+
// Add the constraint with penalization
gmm::mult(gmm::transposed(M4), M4, aux2);
gmm::add (gmm::scaled(aux2, 1.E7), M2);
@@ -832,11 +859,11 @@ namespace getfem {
} else
{ gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
- if (pf1->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof1/Q, Q);
+ if (pf3->target_dim() == 1 && Q > 1) {
+ gmm::sub_slice I(0, ndof3/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
for (size_type i = 0; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof1/Q, Q);
+ gmm::sub_slice I2(i, ndof3/Q, Q);
gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
}
} else
@@ -845,19 +872,17 @@ namespace getfem {
gmm::mult(M2inv, M1, MD);
gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
- // S = (I - inv(M7)*M9)(I - inv(M7)*M8*MD)
- base_matrix MPB(ndof1, ndof1);
- gmm::mult(M7inv, M9, MPB);
- gmm::copy(gmm::identity_matrix(), M9);
- gmm::add(gmm::scaled(MPB, scalar_type(-1)), M9);
-
- base_matrix MPC(ndof1, ndof1), MPD(ndof1, ndof1);
- gmm::mult(M8, MD, MPC);
- gmm::mult(M7inv, MPC, MPD);
- gmm::copy(gmm::identity_matrix(), M7);
- gmm::add(gmm::scaled(MPD, scalar_type(-1)), M7);
-
- gmm::mult(M9, M7, M);
+ // S = (I - inv(M7)*M10)*inv(M7)*(M9 - M8*MD)
+ base_matrix MPB(ndof3, ndof3);
+ gmm::mult(M7inv, M10, MPB);
+ gmm::copy(gmm::identity_matrix(), M10);
+ gmm::add(gmm::scaled(MPB, scalar_type(-1)), M10);
+
+ base_matrix MPC(ndof3, ndof1);
+ gmm::mult(gmm::scaled(M8, scalar_type(-1)), MD, MPC);
+ gmm::add(M9, MPC);
+ gmm::mult(M7inv, MPC, M9);
+ gmm::mult(M10, M9, M);
gmm::clean(M, 1E-13);
}
};
@@ -878,7 +903,7 @@ namespace getfem {
size_type cv, base_matrix &M) const {
// The reconstructed variable "S" is described on mf2 and computed by
// S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)))
- // where P__{\dT} et P_T are L2 projections on the boundary and on the
+ // where P_{\dT} et P_T are L2 projections on the boundary and on the
// interior of T on the corresponding discrete spaces.
// Note that P_{\dT}(v_{dT}) = v_{dT} and P_T(v_T) = v_T and D is
// the reconstructed value on P^{k+1} given by the formula:
@@ -892,21 +917,20 @@ namespace getfem {
// "v_{dT}" is the field on the boundary of T and "n" is the outward
// unit normal.
// The implemented formula is
- // S(v) = v_{dT} - P_{\dT}D(v) - P_{\dT}(v_T) + P_{\dT}(P_T(D(v)))
- // by the mean of the projection matrix from P^{k+1} to the original
space
+ // S(v) = P_{\dT}(v_{dT} - D(v) - P_T(v_T - D(v)) )
+ // by the mean of the projection matrix from P^{k+1} to the target space
// and the projection matrix from interior space to the boundary space.
// As it is built, S(v) is zero on interior dofs.
-
- GMM_ASSERT1(&mf1 == &mf2, "The HHO stabilization transformation is "
- "only defined on the HHO space to itself");
-
+
// Obtaining the fem descriptors
pfem pf1 = mf1.fem_of_element(cv);
short_type degree = pf1->estimated_degree();
bgeot::pgeometric_trans pgt = mf1.linked_mesh().trans_of_convex(cv);
pfem pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be
changed to an
// interior PK method
- pfem pfi = interior_fem_of_hho_method(pf1);
+ pfem pf3 = mf2.fem_of_element(cv);
+ pfem pf1i = interior_fem_of_hho_method(pf1);
+ pfem pf3i = interior_fem_of_hho_method(pf3);
papprox_integration pim
= classical_approx_im(pgt, dim_type(2*degree+2))->approx_method();
@@ -918,11 +942,15 @@ namespace getfem {
= HHO_pgp_pool(pgt, pim->pintegration_points());
pfem_precomp pfp1 = HHO_pfp_pool(pf1, pim->pintegration_points());
pfem_precomp pfp2 = HHO_pfp_pool(pf2, pim->pintegration_points());
- pfem_precomp pfpi = HHO_pfp_pool(pfi, pim->pintegration_points());
+ pfem_precomp pfp3 = HHO_pfp_pool(pf3, pim->pintegration_points());
+ pfem_precomp pfp1i = HHO_pfp_pool(pf1i, pim->pintegration_points());
+ pfem_precomp pfp3i = HHO_pfp_pool(pf3i, pim->pintegration_points());
fem_interpolation_context ctx1(pgp, pfp1, 0, G, cv);
fem_interpolation_context ctx2(pgp, pfp2, 0, G, cv);
- fem_interpolation_context ctxi(pgp, pfpi, 0, G, cv);
+ fem_interpolation_context ctx3(pgp, pfp3, 0, G, cv);
+ fem_interpolation_context ctx1i(pgp, pfp1i, 0, G, cv);
+ fem_interpolation_context ctx3i(pgp, pfp3i, 0, G, cv);
size_type Q = mf1.get_qdim(), N = mf1.linked_mesh().dim();
GMM_ASSERT1(Q == N, "This transformation works only for vector fields "
@@ -932,38 +960,45 @@ namespace getfem {
size_type ndof1 = pf1->nb_dof(cv) * qmult1;
size_type qmult2 = N / pf2->target_dim();
size_type ndof2 = pf2->nb_dof(cv) * qmult2;
- size_type qmulti = N / pfi->target_dim();
- size_type ndofi = pfi->nb_dof(cv) * qmulti;
+ size_type qmult3 = Q / pf3->target_dim();
+ size_type ndof3 = pf3->nb_dof(cv) * qmult3;
+ size_type qmult1i = Q / pf1i->target_dim();
+ size_type ndof1i = pf1i->nb_dof(cv) * qmult1i;
+ size_type qmult3i = Q / pf3i->target_dim();
+ size_type ndof3i = pf3i->nb_dof(cv) * qmult3i;
- base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
- base_matrix tv1p, tv2p, tvi;
+ base_tensor t1, t2, t3, t1i, t3i, tv1, tv2, t1p, t2p;
+ base_matrix tv1p, tv2p, tv3p, tv1i, tv3i;
base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
base_matrix M3(N, ndof1), M4(N, ndof2);
base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
base_matrix M5(N*N, ndof1), M6(N*N, ndof2);
- base_matrix M7(ndof1, ndof1), M7inv(ndof1, ndof1), M8(ndof1, ndof2);
- base_matrix M9(ndof1, ndof1), MD(ndof2, ndof1);
+ base_matrix M7(ndof3, ndof3), M7inv(ndof3, ndof3), M8(ndof3, ndof2);
+ base_matrix M9(ndof3, ndof1), M10(ndof3, ndof3), MD(ndof2, ndof1);
// Integrals on the element : \int_T Sym(Grad(D)).Grad(w) (M2)
// \int_T Sym(Grad(v_T)).Grad(w) (M1)
// \int_T D (M4) and \int_T v_T (M3)
// \int_T Skew(Grad(D)) (M6)
for (size_type ipt = 0; ipt < pim->nb_points_on_convex(); ++ipt) {
- ctx1.set_ii(ipt); ctx2.set_ii(ipt);
+ ctx1.set_ii(ipt); ctx2.set_ii(ipt); ctx3.set_ii(ipt);
scalar_type coeff = pim->coeff(ipt) * ctx1.J();
ctx1.grad_base_value(t1);
- vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), N);
+ vectorize_grad_base_tensor(t1, tv1, ndof1, pf1->target_dim(), Q);
ctx1.base_value(t1p);
- vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), N);
+ vectorize_base_tensor(t1p, tv1p, ndof1, pf1->target_dim(), Q);
ctx2.grad_base_value(t2);
- vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
+ vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
ctx2.base_value(t2p);
- vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), N);
+ vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
+
+ ctx3.base_value(t3);
+ vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
for (size_type i = 0; i < ndof2; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
@@ -995,40 +1030,51 @@ namespace getfem {
for (size_type k = 0; k < N; ++k)
M3(k, i) += coeff * tv1p(i, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
- for (size_type j = 0; j < ndof1; ++j)
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof3; ++j)
for (size_type k = 0; k < N; ++k)
- M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
+ M7(i, j) += coeff * tv3p(i,k) * tv3p(j, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
for (size_type k = 0; k < N; ++k)
- M8(i, j) += coeff * tv1p(i,k) * tv2p(j, k);
+ M8(i, j) += coeff * tv3p(i,k) * tv2p(j, k);
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof1; ++j)
+ for (size_type k = 0; k < Q; ++k)
+ M9(i, j) += coeff * tv3p(i, k) * tv1p(j, k);
}
// Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Grad(w).n) (M1)
// \int_{dT} Skew(n x v_{dT} - v_{dT} x n) (M5)
for (short_type ifc = 0; ifc < pgt->structure()->nb_faces(); ++ifc) {
- ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctxi.set_face_num(ifc);
+ ctx1.set_face_num(ifc); ctx2.set_face_num(ifc); ctx3.set_face_num(ifc);
+ ctx1i.set_face_num(ifc); ctx3i.set_face_num(ifc);
size_type first_ind = pim->ind_first_point_on_face(ifc);
for (size_type ipt = 0; ipt < pim->nb_points_on_face(ifc); ++ipt) {
- ctx1.set_ii(first_ind+ipt);
- ctx2.set_ii(first_ind+ipt);
- ctxi.set_ii(first_ind+ipt);
+ ctx1.set_ii(first_ind+ipt); ctx2.set_ii(first_ind+ipt);
+ ctx3.set_ii(first_ind+ipt); ctx1i.set_ii(first_ind+ipt);
+ ctx3i.set_ii(first_ind+ipt);
scalar_type coeff = pim->coeff(first_ind+ipt) * ctx1.J();
ctx2.grad_base_value(t2);
- vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
+ vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
ctx2.base_value(t2p);
- vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), N);
+ vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
ctx1.base_value(t1);
- vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), N);
+ vectorize_base_tensor(t1, tv1p, ndof1, pf1->target_dim(), Q);
- ctxi.base_value(ti);
- vectorize_base_tensor(ti, tvi, ndofi, pfi->target_dim(), N);
+ ctx1i.base_value(t1i);
+ vectorize_base_tensor(t1i, tv1i, ndof1i, pf1i->target_dim(), Q);
+
+ ctx3i.base_value(t3i);
+ vectorize_base_tensor(t3i, tv3i, ndof3i, pf3i->target_dim(), Q);
+
+ ctx3.base_value(t3);
+ vectorize_base_tensor(t3, tv3p, ndof3, pf3->target_dim(), Q);
gmm::mult(ctx1.B(), pgt->normals()[ifc], un);
@@ -1036,7 +1082,7 @@ namespace getfem {
for (size_type j = 0; j < ndof2; ++j)
for (size_type k1 = 0; k1 < N; ++k1) {
scalar_type b(0), a = coeff *
- (tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
+ (tv1p(i, k1) - (i < ndof1i ? tv1i(i, k1) : 0.));
for (size_type k2 = 0; k2 < N; ++k2)
b += a * 0.5 * (tv2.as_vector()[j+(k1 + k2*N)*ndof2] +
tv2.as_vector()[j+(k2 + k1*N)*ndof2])*
un[k2];
@@ -1049,20 +1095,25 @@ namespace getfem {
M5(k1+k2*N, i) += 0.5 * coeff * (tv1p(i, k1) * un[k2] -
tv1p(i, k2) * un[k1]);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
- for (size_type j = 0; j < ndof1; ++j)
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof3; ++j)
for (size_type k = 0; k < N; ++k)
- M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
+ M7(i, j) += coeff * tv3p(i,k) * tv3p(j, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
for (size_type k = 0; k < N; ++k)
- M8(i, j) += coeff * tv1p(i,k) * tv2p(j, k);
+ M8(i, j) += coeff * tv3p(i,k) * tv2p(j, k);
- for (size_type i = 0; i < ndof1; ++i) // To be optimized
- for (size_type j = 0; j < ndofi; ++j)
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof1; ++j)
+ for (size_type k = 0; k < Q; ++k)
+ M9(i, j) += coeff * tv3p(i, k) * tv1p(j, k);
+
+ for (size_type i = 0; i < ndof3; ++i) // To be optimized
+ for (size_type j = 0; j < ndof3i; ++j)
for (size_type k = 0; k < N; ++k)
- M9(i, j) += coeff * tv1p(i,k) * tvi(j, k);
+ M10(i, j) += coeff * tv3p(i,k) * tv3i(j, k);
}
}
@@ -1078,11 +1129,11 @@ namespace getfem {
gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv);
- if (pf1->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof1/Q, Q);
+ if (pf3->target_dim() == 1 && Q > 1) {
+ gmm::sub_slice I(0, ndof3/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
for (size_type i = 0; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof1/Q, Q);
+ gmm::sub_slice I2(i, ndof3/Q, Q);
gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
}
} else
@@ -1091,19 +1142,17 @@ namespace getfem {
gmm::mult(M2inv, M1, MD);
gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
- // S = (I - inv(M7)*M9)(I - inv(M7)*M8*MD)
- base_matrix MPB(ndof1, ndof1);
- gmm::mult(M7inv, M9, MPB);
- gmm::copy(gmm::identity_matrix(), M9);
- gmm::add(gmm::scaled(MPB, scalar_type(-1)), M9);
-
- base_matrix MPC(ndof1, ndof1), MPD(ndof1, ndof1);
- gmm::mult(M8, MD, MPC);
- gmm::mult(M7inv, MPC, MPD);
- gmm::copy(gmm::identity_matrix(), M7);
- gmm::add(gmm::scaled(MPD, scalar_type(-1)), M7);
-
- gmm::mult(M9, M7, M);
+ // S = (I - inv(M7)*M10)*inv(M7)*(M9 - M8*MD)
+ base_matrix MPB(ndof3, ndof3);
+ gmm::mult(M7inv, M10, MPB);
+ gmm::copy(gmm::identity_matrix(), M10);
+ gmm::add(gmm::scaled(MPB, scalar_type(-1)), M10);
+
+ base_matrix MPC(ndof3, ndof1);
+ gmm::mult(gmm::scaled(M8, scalar_type(-1)), MD, MPC);
+ gmm::add(M9, MPC);
+ gmm::mult(M7inv, MPC, M9);
+ gmm::mult(M10, M9, M);
gmm::clean(M, 1E-13);
}
};
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index fc05256..00d184a 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1666,7 +1666,7 @@ namespace getfem {
/* ******************************************************************** */
// Not tested. To be tested
- struct IPK_SQUARE_ : public fem<base_poly> {
+ struct CIPK_SQUARE_ : public fem<base_poly> {
dim_type k; //
mutable base_matrix K;
base_small_vector norient;
@@ -1678,10 +1678,10 @@ namespace getfem {
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
- IPK_SQUARE_(dim_type nc_);
+ CIPK_SQUARE_(dim_type nc_);
};
- void IPK_SQUARE_::mat_trans(base_matrix &M,
+ void CIPK_SQUARE_::mat_trans(base_matrix &M,
const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
// All the dof of this element are concentrated at the center of the
element
@@ -1730,7 +1730,7 @@ namespace getfem {
}
}
- IPK_SQUARE_::IPK_SQUARE_(dim_type k_) {
+ CIPK_SQUARE_::CIPK_SQUARE_(dim_type k_) {
k = k_;
pgt_stored = 0;
@@ -1762,7 +1762,7 @@ namespace getfem {
}
}
- static pfem IPK_SQUARE(fem_param_list ¶ms,
+ static pfem CIPK_SQUARE(fem_param_list ¶ms,
std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
<< params.size() << " should be 1.");
@@ -1770,7 +1770,7 @@ namespace getfem {
int k = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(k >= 0 && k < 50 && double(k) == params[0].num(),
"Bad parameter");
- pfem p = std::make_shared<IPK_SQUARE_>(dim_type(k));
+ pfem p = std::make_shared<CIPK_SQUARE_>(dim_type(k));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
@@ -3823,10 +3823,9 @@ namespace getfem {
PK_discont_(dim_type nc, short_type k, scalar_type alpha=scalar_type(0))
: PK_fem_(nc, k) {
- if (alpha < 1e-4) // In order to glue dof of two faces in 3D,
- std::fill(dof_types_.begin(), // for alpha > 1e-4. Important.
- dof_types_.end(), lagrange_nonconforming_dof(nc));
-
+ std::fill(dof_types_.begin(),
+ dof_types_.end(), lagrange_nonconforming_dof(nc));
+
if (alpha != scalar_type(0)) {
base_node G =
gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
@@ -3864,6 +3863,139 @@ namespace getfem {
return p;
}
+ /* ******************************************************************** */
+ /* Connectable interior PK */
+ /* ******************************************************************** */
+
+ struct PK_conn_int_ : public PK_fem_ {
+ public :
+
+ PK_conn_int_(dim_type nc, short_type k)
+ : PK_fem_(nc, k) {
+ scalar_type alpha = 0.1;
+
+ std::fill(dof_types_.begin(), dof_types_.end(), lagrange_dof(nc));
+ if (alpha != scalar_type(0)) {
+ base_node G =
+ gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
+ for (size_type i=0; i < cv_node.nb_points(); ++i)
+ cv_node.points()[i] = (1-alpha)*cv_node.points()[i] + alpha*G;
+ for (size_type d = 0; d < nc; ++d) {
+ base_poly S(1,2);
+ S[0] = -alpha * G[d] / (1-alpha);
+ S[1] = 1. / (1-alpha);
+ for (size_type j=0; j < nb_base(0); ++j)
+ base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
+ }
+ }
+ }
+ };
+
+ static pfem conn_int_PK_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
+ << params.size() << " should be 2.");
+ GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
+ "Bad type of parameters");
+ int n = int(::floor(params[0].num() + 0.01));
+ int k = int(::floor(params[1].num() + 0.01));
+ GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150
+ && double(n) == params[0].num()
+ && double(k) == params[1].num(), "Bad parameters");
+ pfem p = std::make_shared<PK_conn_int_>(dim_type(n), short_type(k));
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
+ /* ******************************************************************** */
+ /* Interior PK */
+ /* ******************************************************************** */
+ /* Interior PK element for arbitrary shape */
+ /* (equivalent to DISCONTINUOUS_PK(n,k,0.1) for simplicies */
+
+
+ struct PK_int_ : public PK_fem_ {
+ public :
+
+ PK_int_(dim_type nc, short_type k, bgeot::pconvex_ref cvr_)
+ : PK_fem_(nc, k) {
+ cvr = cvr_;
+
+ scalar_type alpha = 0.1;
+ std::fill(dof_types_.begin(),
+ dof_types_.end(), lagrange_nonconforming_dof(nc));
+
+ if (alpha != scalar_type(0)) {
+ base_node G =
+ gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
+ for (size_type i=0; i < cv_node.nb_points(); ++i)
+ cv_node.points()[i] = (1-alpha)*cv_node.points()[i] + alpha*G;
+ for (size_type d = 0; d < nc; ++d) {
+ base_poly S(1,2);
+ S[0] = -alpha * G[d] / (1-alpha);
+ S[1] = 1. / (1-alpha);
+ for (size_type j=0; j < nb_base(0); ++j)
+ base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
+ }
+ }
+ }
+ };
+
+ static pfem simplex_IPK_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
+ << params.size() << " should be 2.");
+ GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
+ "Bad type of parameters");
+ int n = int(::floor(params[0].num() + 0.01));
+ int k = int(::floor(params[1].num() + 0.01));
+ GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150
+ && double(n) == params[0].num()
+ && double(k) == params[1].num(), "Bad parameters");
+ pfem p = std::make_shared<PK_int_>(dim_type(n), short_type(k),
+ bgeot::simplex_of_reference(dim_type(n)));
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
+ static pfem parallelepiped_IPK_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
+ << params.size() << " should be 2.");
+ GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
+ "Bad type of parameters");
+ int n = int(::floor(params[0].num() + 0.01));
+ int k = int(::floor(params[1].num() + 0.01));
+ GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150
+ && double(n) == params[0].num()
+ && double(k) == params[1].num(), "Bad parameters");
+ pfem p = std::make_shared<PK_int_>(dim_type(n), short_type(k),
+ bgeot::parallelepiped_of_reference(dim_type(n)));
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
+ static pfem prism_IPK_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
+ << params.size() << " should be 2.");
+ GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
+ "Bad type of parameters");
+ int n = int(::floor(params[0].num() + 0.01));
+ int k = int(::floor(params[1].num() + 0.01));
+ GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150
+ && double(n) == params[0].num()
+ && double(k) == params[1].num(), "Bad parameters");
+ pfem p = std::make_shared<PK_int_>(dim_type(n), short_type(k),
+ bgeot::prism_of_reference(dim_type(n)));
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
/* ******************************************************************** */
/* PK element with a bubble base function */
@@ -4052,6 +4184,10 @@ namespace getfem {
add_suffix("PK_DISCONTINUOUS", PK_discontinuous_fem);
add_suffix("PRISM_PK_DISCONTINUOUS", prism_PK_discontinuous_fem);
add_suffix("PK_PRISM_DISCONTINUOUS", prism_PK_discontinuous_fem); // for
backwards compatibility
+ add_suffix("SIMPLEX_IPK", simplex_IPK_fem);
+ add_suffix("PRISM_IPK", prism_IPK_fem);
+ add_suffix("QUAD_IPK", parallelepiped_IPK_fem);
+ add_suffix("SIMPLEX_CIPK", conn_int_PK_fem);
add_suffix("PK_WITH_CUBIC_BUBBLE", PK_with_cubic_bubble);
add_suffix("PRODUCT", product_fem);
add_suffix("P1_NONCONFORMING", P1_nonconforming_fem);
@@ -4068,7 +4204,7 @@ namespace getfem {
add_suffix("PK_FULL_HIERARCHICAL_COMPOSITE",
PK_composite_full_hierarch_fem);
add_suffix("PK_GAUSSLOBATTO1D", PK_GL_fem);
- add_suffix("IPK_QUAD", IPK_SQUARE);
+ add_suffix("QUAD_CIPK", CIPK_SQUARE);
add_suffix("Q2_INCOMPLETE", Q2_incomplete_fem);
add_suffix("Q2_INCOMPLETE_DISCONTINUOUS",
Q2_incomplete_discontinuous_fem);
add_suffix("HCT_TRIANGLE", HCT_triangle_fem);
- [Getfem-commits] [getfem-commits] master updated (860f866 -> e920df7), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject),
Yves Renard <=
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27