[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:03 -0400 (EDT) |
branch: master
commit 955324bbd996e3889c14d9aeb5fafd548307ddde
Author: Yves Renard <address@hidden>
Date: Mon Aug 19 16:58:33 2019 +0200
Adding a test of HHO methods on a Poisson problem and various fixes
---
interface/src/gf_mesh_get.cc | 32 +++
interface/src/gf_model_set.cc | 63 ++++++
interface/tests/python/Makefile.am | 6 +-
interface/tests/python/demo_laplacian.py | 4 +-
.../{demo_laplacian.py => demo_laplacian_HHO.py} | 68 ++++++-
src/bgeot_geometric_trans.cc | 2 +-
src/bgeot_poly_composite.cc | 76 ++++----
src/getfem/bgeot_poly.h | 2 +-
src/getfem/bgeot_poly_composite.h | 29 ++-
src/getfem/getfem_HHO.h | 16 +-
src/getfem/getfem_fem.h | 5 +-
src/getfem/getfem_mesh.h | 7 +
src/getfem_HHO.cc | 216 +++++++++++----------
src/getfem_fem.cc | 2 +-
src/getfem_fem_composite.cc | 8 +-
src/getfem_generic_assembly_compile_and_exec.cc | 16 +-
src/getfem_generic_assembly_semantic.cc | 12 +-
src/getfem_mesh.cc | 17 ++
18 files changed, 393 insertions(+), 188 deletions(-)
diff --git a/interface/src/gf_mesh_get.cc b/interface/src/gf_mesh_get.cc
index 20dc2da..077861d 100644
--- a/interface/src/gf_mesh_get.cc
+++ b/interface/src/gf_mesh_get.cc
@@ -244,6 +244,28 @@ outer_faces(const getfem::mesh &m, mexargs_in &in,
mexargs_out &out, const std::
}
+static void
+all_faces(const getfem::mesh &m, mexargs_in &in, mexargs_out &out) {
+ dal::bit_vector cvlst;
+
+ if (in.remaining()) cvlst = in.pop().to_bit_vector(&m.convex_index());
+ else cvlst = m.convex_index();
+ getfem::mesh_region mr;
+ for (dal::bv_visitor ic(cvlst); !ic.finished(); ++ic) mr.add(ic);
+ getfem::mesh_region mrr = all_faces_of_mesh(m, mr);
+
+ unsigned fcnt = 0;
+ for (getfem::mr_visitor i(mrr); !i.finished(); ++i) ++fcnt;
+ iarray w = out.pop().create_iarray(2, fcnt);
+ fcnt = 0;
+ for (getfem::mr_visitor i(mrr); !i.finished(); ++i) {
+ w(0,fcnt) = int(i.cv()+config::base_index());
+ w(1,fcnt) = int(short_type(i.f()+config::base_index()));
+ fcnt++;
+ }
+}
+
+
static void
inner_faces(const getfem::mesh &m, mexargs_in &in, mexargs_out &out) {
@@ -789,6 +811,16 @@ void gf_mesh_get(getfemint::mexargs_in& m_in,
check_empty_mesh(pmesh);
inner_faces(*pmesh, in, out);
);
+
+ /*@GET CVFIDs = ('all faces'[, CVIDs])
+ Return the set of faces of the in CVIDs (in all the mesh if CVIDs is
+ omitted). Note that the face shared by two neighbour elements will be
+ represented twice. @*/
+ sub_command
+ ("all faces", 0, 1, 0, 1,
+ check_empty_mesh(pmesh);
+ all_faces(*pmesh, in, out);
+ );
/*@GET CVFIDs = ('outer faces with direction', @vec v, @scalar angle [,
CVIDs])
Return the set of faces not shared by two convexes and with a mean outward
vector lying within an angle `angle` (in radians) from vector `v`.
diff --git a/interface/src/gf_model_set.cc b/interface/src/gf_model_set.cc
index 5846e78..8f1ad05 100644
--- a/interface/src/gf_model_set.cc
+++ b/interface/src/gf_model_set.cc
@@ -29,6 +29,7 @@
#include <getfem/getfem_contact_and_friction_large_sliding.h>
#include <getfem/getfem_nonlinear_elasticity.h>
#include <getfem/getfem_plasticity.h>
+#include <getfem/getfem_HHO.h>
#include <getfem/getfem_fourth_order.h>
#include <getfem/getfem_linearized_plates.h>
#include <getfemint_misc.h>
@@ -443,6 +444,68 @@ void gf_model_set(getfemint::mexargs_in& m_in,
add_P0_projection(*md, transname);
);
+ /*@SET ('add HHO reconstructed gradient', @str transname)
+ Add to the model the elementary transformation corresponding to the
+ reconstruction of a gradient for HHO methods.
+ The name is the name given to the elementary transformation. @*/
+ sub_command
+ ("add HHO reconstructed gradient", 1, 1, 0, 0,
+ std::string transname = in.pop().to_string();
+ add_HHO_reconstructed_gradient(*md, transname);
+ );
+
+ /*@SET ('add HHO reconstructed symmetrized gradient', @str transname)
+ Add to the model the elementary transformation corresponding to the
+ reconstruction of a symmetrized gradient for HHO methods.
+ The name is the name given to the elementary transformation. @*/
+ sub_command
+ ("add HHO reconstructed symmetrized gradient", 1, 1, 0, 0,
+ std::string transname = in.pop().to_string();
+ add_HHO_reconstructed_symmetrized_gradient(*md, transname);
+ );
+
+ /*@SET ('add HHO reconstructed value', @str transname)
+ Add to the model the elementary transformation corresponding to the
+ reconstruction of the variable for HHO methods.
+ The name is the name given to the elementary transformation. @*/
+ sub_command
+ ("add HHO reconstructed value", 1, 1, 0, 0,
+ std::string transname = in.pop().to_string();
+ add_HHO_reconstructed_value(*md, transname);
+ );
+
+ /*@SET ('add HHO reconstructed symmetrized value', @str transname)
+ Add to the model the elementary transformation corresponding to the
+ reconstruction of the variable for HHO methods using a symmetrized
+ gradient.
+ The name is the name given to the elementary transformation. @*/
+ sub_command
+ ("add HHO reconstructed symmetrized value", 1, 1, 0, 0,
+ std::string transname = in.pop().to_string();
+ add_HHO_reconstructed_symmetrized_value(*md, transname);
+ );
+
+ /*@SET ('add HHO stabilization', @str transname)
+ Add to the model the elementary transformation corresponding to the
+ HHO stabilization operator.
+ The name is the name given to the elementary transformation. @*/
+ sub_command
+ ("add HHO stabilization", 1, 1, 0, 0,
+ std::string transname = in.pop().to_string();
+ add_HHO_stabilization(*md, transname);
+ );
+
+ /*@SET ('add HHO symmetrized stabilization', @str transname)
+ Add to the model the elementary transformation corresponding to the
+ HHO stabilization operator using a symmetrized gradient.
+ The name is the name given to the elementary transformation. @*/
+ sub_command
+ ("add HHO symmetrized stabilization", 1, 1, 0, 0,
+ std::string transname = in.pop().to_string();
+ add_HHO_symmetrized_stabilization(*md, transname);
+ );
+
+
/*@SET ('add interpolate transformation from expression', @str transname,
@tmesh source_mesh, @tmesh target_mesh, @str expr)
Add a transformation to the model from mesh `source_mesh` to mesh
`target_mesh` given by the expression `expr` which corresponds to a
diff --git a/interface/tests/python/Makefile.am
b/interface/tests/python/Makefile.am
index cd6ff0a..113a191 100644
--- a/interface/tests/python/Makefile.am
+++ b/interface/tests/python/Makefile.am
@@ -1,4 +1,4 @@
-# Copyright (C) 1999-2018 Yves Renard
+# Copyright (C) 1999-2019 Yves Renard
#
# This file is a part of GetFEM++
#
@@ -32,6 +32,7 @@ EXTRA_DIST= \
demo_crack.py \
demo_fictitious_domains.py \
demo_laplacian.py \
+ demo_laplacian_HHO.py \
demo_laplacian_pyramid.py \
demo_laplacian_DG.py \
demo_laplacian_aposteriori.py \
@@ -39,7 +40,7 @@ EXTRA_DIST= \
demo_phase_field.py \
demo_plasticity.py \
demo_plate.py \
- demo_unit_disk.py \
+ demo_unit_disk.py \
demo_static_contact.py \
demo_dynamic_contact_1D.py \
demo_Mindlin_Reissner_plate.py \
@@ -68,6 +69,7 @@ TESTS = \
demo_wave.py \
demo_wave_equation.py \
demo_laplacian.py \
+ demo_laplacian_HHO.py \
$(optpy)
AM_TESTS_ENVIRONMENT = \
diff --git a/interface/tests/python/demo_laplacian.py
b/interface/tests/python/demo_laplacian.py
index b2da830..6e803b6 100644
--- a/interface/tests/python/demo_laplacian.py
+++ b/interface/tests/python/demo_laplacian.py
@@ -2,7 +2,7 @@
# -*- coding: utf-8 -*-
# Python GetFEM++ interface
#
-# Copyright (C) 2004-2017 Yves Renard, Julien Pommier.
+# Copyright (C) 2004-2019 Yves Renard, Julien Pommier.
#
# This file is a part of GetFEM++
#
@@ -43,7 +43,7 @@ m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
# Create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
mfu = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
-# assign the P2 fem to all convexes of the both MeshFem
+# assign the P2 fem to all elements of the both MeshFem
mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
mfrhs.set_fem(gf.Fem('FEM_PK(2,2)'))
diff --git a/interface/tests/python/demo_laplacian.py
b/interface/tests/python/demo_laplacian_HHO.py
similarity index 69%
copy from interface/tests/python/demo_laplacian.py
copy to interface/tests/python/demo_laplacian_HHO.py
index b2da830..f091623 100644
--- a/interface/tests/python/demo_laplacian.py
+++ b/interface/tests/python/demo_laplacian_HHO.py
@@ -2,7 +2,7 @@
# -*- coding: utf-8 -*-
# Python GetFEM++ interface
#
-# Copyright (C) 2004-2017 Yves Renard, Julien Pommier.
+# Copyright (C) 2019-2019 Yves Renard.
#
# This file is a part of GetFEM++
#
@@ -19,7 +19,7 @@
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
-""" 2D Poisson problem test.
+""" 2D 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++.
@@ -35,16 +35,27 @@ NX = 100 # Mesh parameter.
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
+using_HHO = False # 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();
-# Create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
+# Meshfems
mfu = gf.MeshFem(m, 1)
+mfgu = gf.MeshFem(m, N)
+mfur = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
-# assign the P2 fem to all convexes of the both MeshFem
-mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
+
+if (using_HHO):
+ mfu.set_fem(gf.Fem('FEM_HHO(FEM_PK(2,2),FEM_PK(1,2))'))
+ mfgu.set_fem(gf.Fem('FEM_HHO(FEM_PK(2,2),FEM_PK(1,2))'))
+else:
+ mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
+ mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
+
+mfur.set_fem(gf.Fem('FEM_PK(2,3)'))
mfrhs.set_fem(gf.Fem('FEM_PK(2,2)'))
# Integration method used
@@ -67,6 +78,11 @@ m.set_region(DIRICHLET_BOUNDARY_NUM1, fleft)
m.set_region(DIRICHLET_BOUNDARY_NUM2, ftop)
m.set_region(NEUMANN_BOUNDARY_NUM, fneum)
+# Faces for stabilization term
+all_faces = m.all_faces()
+ALL_FACES = 4
+m.set_region(ALL_FACES, all_faces)
+
# Interpolate the exact solution (Assuming mfu is a Lagrange fem)
Ue = mfu.eval('y*(y-1)*x*(x-1)+x*x*x*x*x')
@@ -77,11 +93,35 @@ F2 = mfrhs.eval('[y*(y-1)*(2*x-1) + 5*x*x*x*x,
x*(x-1)*(2*y-1)]')
# Model
md = gf.Model('real')
+
+
# Main unknown
md.add_fem_variable('u', mfu)
+md.add_fem_data('Gu', mfgu)
+md.add_fem_data('ur', mfur)
+
+# Needed reconstuction and stabilization operators
+if (using_HHO):
+ md.add_HHO_reconstructed_gradient('HHO_Grad');
+ md.add_HHO_reconstructed_value('HHO_Val');
+ md.add_HHO_stabilization('HHO_Stab');
+ md.add_macro('HHO_Val_u', 'Elementary_transformation(u, HHO_Val)')
+ md.add_macro('HHO_Grad_u', 'Elementary_transformation(u, HHO_Grad, Gu)')
+ md.add_macro('HHO_Grad_Test_u',
+ 'Elementary_transformation(Test_u, HHO_Grad, Gu)')
+ md.add_macro('HHO_Stab_u', 'Elementary_transformation(u, HHO_Stab)')
+ md.add_macro('HHO_Stab_Test_u',
+ 'Elementary_transformation(Test_u, HHO_Stab, ur)')
+
# Laplacian term on u
-md.add_Laplacian_brick(mim, 'u')
+if (using_HHO):
+ # Laplacian term
+ md.add_linear_term(mim, 'HHO_Grad_u.HHO_Grad_Test_u')
+ # Stabilization term
+ md.add_linear_term(mim, 'HHO_Stab_u.HHO_Stab_Test_u', ALL_FACES)
+else:
+ md.add_Laplacian_brick(mim, 'u')
# Volumic source term
md.add_initialized_fem_data('VolumicData', mfrhs, F1)
@@ -124,11 +164,23 @@ md.solve()
# Main unknown
U = md.variable('u')
-L2error = gf.compute(mfu, U-Ue, 'L2 norm', mim)
-H1error = gf.compute(mfu, U-Ue, 'H1 norm', mim)
+if (using_HHO):
+ L2error = gf.asm('generic', mim, 0, 'sqr(HHO_Val_u-Ue)',
+ -1, md, 'Ue', 0, mfu, Ue)
+ H1error = gf.asm('generic', mim, 0, 'Norm_sqr(Grad_u-Grad_Ue)',
+ -1, md, 'Ue', 0, mfu, Ue)
+ H1error = np.sqrt(L2error + H1error)
+ L2error = np.sqrt(L2error)
+else :
+ L2error = gf.compute(mfu, U-Ue, 'L2 norm', mim)
+ H1error = gf.compute(mfu, U-Ue, 'H1 norm', mim)
print('Error in L2 norm : ', L2error)
print('Error in H1 norm : ', H1error)
+
+
+# Calculer l'erreur sur la reconstruction, aussi.
+
# Export data
mfu.export_to_pos('laplacian.pos', Ue,'Exact solution',
U,'Computed solution')
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index b33a92d..363ade1 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -1390,7 +1390,7 @@ namespace bgeot {
}
void delete_geotrans_precomp(pgeotrans_precomp pgp)
- { dal::del_stored_object(pgp); }
+ { dal::del_stored_object(pgp, true); }
} /* end of namespace bgeot. */
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index e37b33c..716d71f 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -81,21 +81,16 @@ namespace bgeot {
pgeometric_trans pgt = m.trans_of_convex(cv);
size_type N = pgt->structure()->dim();
size_type P = m.dim();
- GMM_ASSERT1(pgt->is_linear() && N == P, "Bad geometric transformation");
+ GMM_ASSERT1(pgt->is_linear(), "Bad geometric transformation");
base_matrix G(P, pgt->nb_points());
- base_matrix pc(pgt->nb_points() , N);
- base_matrix B0(N, P);
-
+ base_node X(N);
+
m.points_of_convex(cv, G);
-
- base_node x(N); gmm::clear(x);
- pgt->poly_vector_grad(x, pc);
-
- gmm::mult(gmm::transposed(pc), gmm::transposed(G), B0);
- det[cv] = gmm::lu_inverse(B0);
- gtrans[cv] = B0;
-
+ bgeot::geotrans_interpolation_context ctx(pgt, X, G);
+ gtrans[cv] = ctx.B();
+ det[cv] = ctx.J();
+
auto points_of_convex = m.points_of_convex(cv);
orgs[cv] = points_of_convex[0];
bounding_box(min, max, points_of_convex);
@@ -105,29 +100,35 @@ namespace bgeot {
box_tree.build_tree();
}
- DAL_TRIPLE_KEY(base_poly_key, short_type, short_type,
std::vector<opt_long_scalar_type>);
+ DAL_TRIPLE_KEY(base_poly_key, short_type, short_type,
+ std::vector<opt_long_scalar_type>);
polynomial_composite::polynomial_composite(const mesh_precomposite &m,
bool lc, bool ff)
: mp(&m), local_coordinate(lc), faces_first(ff),
- default_poly(mp->dim(), 0) {}
+ default_polys(mp->dim()+1) {
+ for (dim_type i = 0; i <= mp->dim(); ++i)
+ default_polys[i] = base_poly(i, 0.);
+ }
void polynomial_composite::derivative(short_type k) {
if (local_coordinate) {
- dim_type N = mp->dim();
- base_poly P(N, 0), Q;
- base_vector e(N), f(N);
+ dim_type P = mp->dim();
+ base_vector e(P), f(P); e[k] = 1.0;
for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
- gmm::clear(e); e[k] = 1.0;
+ dim_type N = dim_type(gmm::mat_ncols(mp->gtrans[ic]));
+ f.resize(N);
gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
- P.clear();
- auto &poly = poly_of_subelt(ic);
- for (dim_type n = 0; n < N; ++n) {
- Q = poly;
- Q.derivative(n);
- P += Q * f[n];
+ base_poly Der(N, 0), Q;
+ if (polytab.find(ic) != polytab.end()) {
+ auto &poly = poly_of_subelt(ic);
+ for (dim_type n = 0; n < N; ++n) {
+ Q = poly;
+ Q.derivative(n);
+ Der += Q * f[n];
+ }
+ if (Der.is_zero()) polytab.erase(ic); else
set_poly_of_subelt(ic,Der);
}
- if (polytab.find(ic) != polytab.end()) set_poly_of_subelt(ic, P);
}
}
else
@@ -137,22 +138,29 @@ namespace bgeot {
if (polytab.find(ic) != polytab.end()) set_poly_of_subelt(ic, poly);
}
}
-
- void polynomial_composite::set_poly_of_subelt(size_type l, const base_poly
&poly) {
- auto poly_key = std::make_shared<base_poly_key>(poly.degree(), poly.dim(),
poly);
- auto o = dal::search_stored_object(poly_key);
+
+ void polynomial_composite::set_poly_of_subelt(size_type l,
+ const base_poly &poly) {
+ auto poly_key =
+ std::make_shared<base_poly_key>(poly.degree(), poly.dim(), poly);
+ pstored_base_poly o(std::dynamic_pointer_cast<const stored_base_poly>
+ (dal::search_stored_object(poly_key)));
if (!o) {
- o = std::make_shared<base_poly>(poly);
+ o = std::make_shared<stored_base_poly>(poly);
dal::add_stored_object(poly_key, o);
}
- polytab[l] = poly_key;
+ polytab[l] = o;
}
const base_poly &polynomial_composite::poly_of_subelt(size_type l) const {
auto it = polytab.find(l);
- if (it == polytab.end()) return default_poly;
- return dynamic_cast<const base_poly &>(
- *dal::search_stored_object_on_all_threads(it->second));
+ if (it == polytab.end()) {
+ if (local_coordinate)
+ return default_polys[gmm::mat_ncols(mp->gtrans[l])];
+ else
+ return default_polys[mp->dim()];
+ }
+ return *(it->second);
}
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index ada51fe..1c8dfde 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -177,7 +177,7 @@ namespace bgeot
* The resulting polynomials have a smaller degree.
*
*/
- template<typename T> class polynomial : public std::vector<T>, public
dal::static_stored_object {
+ template<typename T> class polynomial : public std::vector<T> {
protected :
short_type n, d;
diff --git a/src/getfem/bgeot_poly_composite.h
b/src/getfem/bgeot_poly_composite.h
index e704282..744a8c7 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -94,16 +94,22 @@ namespace bgeot {
typedef const mesh_precomposite *pmesh_precomposite;
+ struct stored_base_poly : base_poly, public dal::static_stored_object {
+ stored_base_poly(const base_poly &p) : base_poly(p) {}
+ };
+ typedef std::shared_ptr<const stored_base_poly> pstored_base_poly;
+
+
class polynomial_composite {
protected :
const mesh_precomposite *mp;
- std::map<size_type, dal::pstatic_stored_object_key> polytab;
+ std::map<size_type, pstored_base_poly> polytab;
bool local_coordinate; // Local coordinates on each sub-element for
// polynomials or global coordinates ?
bool faces_first; // If true try to evaluate on faces before on the
// interior, usefull for HHO elements.
- base_poly default_poly;
+ std::vector<base_poly> default_polys;
public :
@@ -147,14 +153,14 @@ namespace bgeot {
if (l != size_type(-1)) {
if (!local_coordinate) return poly_of_subelt(l).eval(it);
- base_node p(mp->dim());
- std::copy(it, it + mp->dim(), p.begin());
+ base_node p(gmm::mat_ncols(mp->gtrans[l]));
+ // std::copy(it, it + mp->dim(), p.begin());
mult_diff_transposed(mp->gtrans[l], it, mp->orgs[l], p);
return poly_of_subelt(l).eval(p.begin());
}
auto dim = mp->dim();
- base_node p0(dim), p0_stored;
+ base_node p0(dim), p1_stored, p1;
size_type cv_stored(-1);
std::copy(it, it + mp->dim(), p0.begin());
@@ -182,12 +188,15 @@ namespace bgeot {
}
for (auto cv : mp->box_to_convexes_map.at(pmax_box->id)) {
- mult_diff_transposed(mp->gtrans[cv], it, mp->orgs[cv], p0);
- if (mp->trans_of_convex(cv)->convex_ref()->is_in(p0) < 1E-10) {
+ cout << "cv = " << cv << endl;
+ p1.resize(gmm::mat_ncols(mp->gtrans[cv]));
+ mult_diff_transposed(mp->gtrans[cv], it, mp->orgs[cv], p1);
+ cout << "p1 = " << p1 << endl;
+ if (mp->trans_of_convex(cv)->convex_ref()->is_in(p1) < 1E-10) {
if (!faces_first || mp->trans_of_convex(cv)->structure()->dim() <
dim)
return to_scalar(poly_of_subelt(cv).eval(local_coordinate
- ? p0.begin() : it));
- p0_stored = p0; cv_stored = cv;
+ ? p1.begin() : it));
+ p1_stored = p1; cv_stored = cv;
}
}
@@ -196,7 +205,7 @@ namespace bgeot {
}
if (cv_stored != size_type(-1))
return to_scalar(poly_of_subelt(cv_stored).eval(local_coordinate
- ? p0_stored.begin():
it));
+ ? p1_stored.begin():
it));
GMM_ASSERT1(false, "Element not found in composite polynomial: "
<< base_node(*it, *it + mp->dim()));
}
diff --git a/src/getfem/getfem_HHO.h b/src/getfem/getfem_HHO.h
index 6eefcbd..d580f5c 100644
--- a/src/getfem/getfem_HHO.h
+++ b/src/getfem/getfem_HHO.h
@@ -42,38 +42,38 @@
namespace getfem {
- /** Add the elementary transformation corresponding to the reconstruction
- of a gradient for HHO methods to the model.
+ /** Add to the model the elementary transformation corresponding to the
+ reconstruction of a gradient for HHO methods.
The name is the name given to the elementary transformation.
*/
void add_HHO_reconstructed_gradient(model &md, std::string name);
- /** Add the elementary transformation corresponding to the reconstruction
- of a symmetrized gradient for HHO methods to the model.
+ /** Add to the model the elementary transformation corresponding to
+ the reconstruction of a symmetrized gradient for HHO methods.
The name is the name given to the elementary transformation.
*/
void add_HHO_reconstructed_symmetrized_gradient(model &md, std::string name);
- /** Add the elementary transformation to the model corresponding to the
+ /** Add to the model the elementary transformation corresponding to the
reconstruction of the variable.
The name is the name given to the elementary transformation.
*/
void add_HHO_reconstructed_value(model &md, std::string name);
- /** Add the elementary transformation to the model corresponding to the
+ /** Add to the model the elementary transformation corresponding to the
reconstruction of the variable using a symmetrized gradient.
The name is the name given to the elementary transformation.
*/
void add_HHO_reconstructed_symmetrized_value(model &md, std::string name);
- /** Add the elementary transformation to the model corresponding to the
+ /** Add to the model the elementary transformation corresponding to the
HHO stabilization operator.
The name is the name given to the elementary transformation.
*/
void add_HHO_stabilization(model &md, std::string name);
- /** Add the elementary transformation to the model corresponding to the
+ /** Add to the model the elementary transformation corresponding to the
HHO stabilization operator using a symmetrized gradient.
The name is the name given to the elementary transformation.
*/
diff --git a/src/getfem/getfem_fem.h b/src/getfem/getfem_fem.h
index 97317c5..295a588 100644
--- a/src/getfem/getfem_fem.h
+++ b/src/getfem/getfem_fem.h
@@ -988,7 +988,10 @@ namespace getfem {
}
}
- /* Specific function for a HHO method to obtain the method in the interior */
+ /**
+ Specific function for a HHO method to obtain the method in the interior.
+ If the method is not of composite type, return the argument.
+ */
pfem interior_fem_of_hho_method(pfem hho_method);
diff --git a/src/getfem/getfem_mesh.h b/src/getfem/getfem_mesh.h
index 1c46be5..46b7edd 100644
--- a/src/getfem/getfem_mesh.h
+++ b/src/getfem/getfem_mesh.h
@@ -654,6 +654,13 @@ namespace getfem {
inner_faces_of_mesh(const mesh &m,
const mesh_region &mr = mesh_region::all_convexes());
+ /** Select all the faces of the given mesh region. The faces are represented*
+ twice if they are shared by two neighbour elements.
+ */
+ mesh_region APIDECL
+ all_faces_of_mesh(const mesh &m,
+ const mesh_region &mr = mesh_region::all_convexes());
+
/** Select in the region mr the faces of the mesh m with their unit
outward vector having a maximal angle "angle" with the vector V.
*/
diff --git a/src/getfem_HHO.cc b/src/getfem_HHO.cc
index 274bd57..2274741 100644
--- a/src/getfem_HHO.cc
+++ b/src/getfem_HHO.cc
@@ -28,6 +28,13 @@ namespace getfem {
THREAD_SAFE_STATIC bgeot::geotrans_precomp_pool HHO_pgp_pool;
THREAD_SAFE_STATIC fem_precomp_pool HHO_pfp_pool;
+ // To be optimized:
+ // - The fact that (when pf2->target_dim() = 1) the
+ // 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
+
+
class _HHO_reconstructed_gradient
: public virtual_elementary_transformation {
@@ -43,12 +50,6 @@ namespace getfem {
// inside the element whose gradient is to be reconstructed,
// "v_{dT}" is the field on the boundary of T and "n" is the outward
// unit normal.
-
- // To be optimized:
- // - The fact that (when pf2->target_dim() = 1) the
- // 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
// Obtaining the fem descriptors
pfem pf1 = mf1.fem_of_element(cv);
@@ -86,7 +87,7 @@ namespace getfem {
base_tensor t1, t2, ti, tv1;
base_matrix tv2, tv1p, tvi;
- base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2);
+ base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
base_matrix aux2(ndof2, ndof2);
// Integrals on the element : \int_T G.w (M2) and \int_T Grad(v_T).w
(M1)
@@ -140,19 +141,21 @@ namespace getfem {
M1(j, i) += b;
}
}
-
}
+
if (pf2->target_dim() == 1) {
- gmm::sub_slice I(0, ndof2, N*Q);
+ gmm::sub_slice I(0, ndof2/(N*Q), N*Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- for (size_type i = 1; i < N*Q; ++i) {
- gmm::sub_slice I2(i, ndof2, N*Q);
- gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2, I2, I2));
+ for (size_type i = 0; i < N*Q; ++i) {
+ gmm::sub_slice I2(i, ndof2/(N*Q), N*Q);
+ gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
}
- } else
- gmm::lu_inverse(M2);
+ } else { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
- gmm::mult(M2, M1, M);
+
+ gmm::mult(M2inv, M1, M);
+ gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
+ // cout << "M = " << M << endl;
}
};
@@ -218,7 +221,7 @@ namespace getfem {
base_tensor t1, t2, ti, tv1;
base_matrix tv2, tv1p, tvi;
- base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2);
+ base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);
base_matrix aux2(ndof2, ndof2);
// Integrals on the element : \int_T G:w (M2)
@@ -277,16 +280,19 @@ namespace getfem {
}
}
+
if (pf2->target_dim() == 1) {
- gmm::sub_slice I(0, ndof2, N*Q);
+ gmm::sub_slice I(0, ndof2/(N*Q), N*Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- for (size_type i = 1; i < N*Q; ++i) {
- gmm::sub_slice I2(i, ndof2, N*Q);
- gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2, I2, I2));
+ for (size_type i = 0; i < N*Q; ++i) {
+ gmm::sub_slice I2(i, ndof2/(N*Q), N*Q);
+ gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
}
} else
- gmm::lu_inverse(M2);
- gmm::mult(M2, M1, M);
+ { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
+
+ gmm::mult(M2inv, M1, M);
+ gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
}
};
@@ -352,7 +358,7 @@ namespace getfem {
base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
base_matrix tv1p, tv2p, tvi;
- base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2);
+ 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);
@@ -373,7 +379,7 @@ namespace getfem {
vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
ctx2.base_value(t2p);
- vectorize_base_tensor(t2p, tv2p, ndof1, pf1->target_dim(), Q);
+ vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
for (size_type i = 0; i < ndof2; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
@@ -433,20 +439,22 @@ namespace getfem {
// Add the constraint with penalization
gmm::mult(gmm::transposed(M4), M4, aux2);
- gmm::add (aux2, M2);
+ gmm::add (gmm::scaled(aux2, 1.E7), M2);
gmm::mult(gmm::transposed(M4), M3, aux1);
- gmm::add (aux1, M1);
+ gmm::add (gmm::scaled(aux1, 1.E7), M1);
if (pf2->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof2, Q);
+ gmm::sub_slice I(0, ndof2/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- for (size_type i = 1; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof2, Q);
- gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2, I2, I2));
+ for (size_type i = 0; i < Q; ++i) {
+ gmm::sub_slice I2(i, ndof2/Q, Q);
+ gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
}
} else
- gmm::lu_inverse(M2);
- gmm::mult(M2, M1, M);
+ { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
+
+ gmm::mult(M2inv, M1, M);
+ gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
}
};
@@ -514,7 +522,7 @@ namespace getfem {
base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
base_matrix tv1p, tv2p, tvi;
- base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2);
+ base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2), M2inv(ndof2, ndof2);;
base_matrix M3(N, ndof1), M4(N, ndof2);
base_matrix M5(N*N, ndof1), M6(N*N, ndof2);
base_matrix aux1(ndof2, ndof1), aux2(ndof2, ndof2);
@@ -537,7 +545,7 @@ namespace getfem {
vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), N);
ctx2.base_value(t2p);
- vectorize_base_tensor(t2p, tv2p, ndof1, pf1->target_dim(), N);
+ vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), N);
for (size_type i = 0; i < ndof2; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
@@ -614,24 +622,26 @@ namespace getfem {
// Add the constraint with penalization
gmm::mult(gmm::transposed(M4), M4, aux2);
- gmm::add (aux2, M2);
+ gmm::add (gmm::scaled(aux2, 1.E7), M2);
gmm::mult(gmm::transposed(M6), M6, aux2);
- gmm::add (aux2, M2);
+ gmm::add (gmm::scaled(aux2, 1.E7), M2);
gmm::mult(gmm::transposed(M4), M3, aux1);
- gmm::add (aux1, M1);
+ gmm::add (gmm::scaled(aux1, 1.E7), M1);
gmm::mult(gmm::transposed(M6), M5, aux1);
- gmm::add (aux1, M1);
+ gmm::add (gmm::scaled(aux1, 1.E7), M1);
if (pf2->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof2, Q);
+ gmm::sub_slice I(0, ndof2/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- for (size_type i = 1; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof2, Q);
- gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2, I2, I2));
+ for (size_type i = 0; i < Q; ++i) {
+ gmm::sub_slice I2(i, ndof2/Q, Q);
+ gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
}
} else
- gmm::lu_inverse(M2);
- gmm::mult(M2, M1, M);
+ { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
+
+ gmm::mult(M2inv, M1, M);
+ gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
}
};
@@ -671,13 +681,13 @@ namespace getfem {
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 pf2 = classical_fem(pgt, short_type(degree + 1)); // Should be
+ // changed for an interior PK method
pfem pfi = interior_fem_of_hho_method(pf1);
papprox_integration pim
@@ -708,10 +718,10 @@ namespace getfem {
base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
base_matrix tv1p, tv2p, tvi;
- base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2);
+ 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), M8(ndof1, ndof2);
+ base_matrix M7(ndof1, ndof1), M7inv(ndof1, ndof1), M8(ndof1, ndof2);
base_matrix M9(ndof1, ndof1), MD(ndof2, ndof1);
// Integrals on the element : \int_T Grad(D).Grad(w) (M2)
@@ -731,7 +741,7 @@ namespace getfem {
vectorize_grad_base_tensor(t2, tv2, ndof2, pf2->target_dim(), Q);
ctx2.base_value(t2p);
- vectorize_base_tensor(t2p, tv2p, ndof1, pf1->target_dim(), Q);
+ vectorize_base_tensor(t2p, tv2p, ndof2, pf2->target_dim(), Q);
for (size_type i = 0; i < ndof2; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
@@ -755,13 +765,13 @@ namespace getfem {
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof1; ++j)
- for (size_type k = 0; k < Q*N; ++k)
- M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
-
+ for (size_type k = 0; k < Q; ++k)
+ M7(i, j) += coeff * tv1p(i, k) * tv1p(j, k);
+
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
- for (size_type k = 0; k < Q*N; ++k)
- M8(i, j) += coeff * tv1p(i,k) * tv2p(j, k);
+ for (size_type k = 0; k < Q; ++k)
+ M8(i, j) += coeff * tv1p(i, k) * tv2p(j, k);
}
@@ -801,64 +811,64 @@ namespace getfem {
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof1; ++j)
- for (size_type k = 0; k < Q*N; ++k)
+ for (size_type k = 0; k < Q; ++k)
M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
- for (size_type k = 0; k < Q*N; ++k)
+ for (size_type k = 0; k < Q; ++k)
M8(i, j) += coeff * tv1p(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 k = 0; k < Q*N; ++k)
- M9(i, j) += coeff * tv1p(i,k) * tvi(j, k);
-
+ for (size_type k = 0; k < Q; ++k)
+ M9(i, j) += coeff * tv1p(i,k) * tvi(j, k);
}
-
}
// Add the constraint with penalization
gmm::mult(gmm::transposed(M4), M4, aux2);
- gmm::add (aux2, M2);
+ gmm::add (gmm::scaled(aux2, 1.E7), M2);
gmm::mult(gmm::transposed(M4), M3, aux1);
- gmm::add (aux1, M1);
+ gmm::add (gmm::scaled(aux1, 1.E7), M1);
if (pf2->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof2, Q);
+ gmm::sub_slice I(0, ndof2/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- for (size_type i = 1; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof2, Q);
- gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2, I2, I2));
+ for (size_type i = 0; i < Q; ++i) {
+ gmm::sub_slice I2(i, ndof2/Q, Q);
+ gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
}
} else
- gmm::lu_inverse(M2);
+ { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
if (pf1->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof1, Q);
+ gmm::sub_slice I(0, ndof1/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
- for (size_type i = 1; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof1, Q);
- gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7, I2, I2));
+ for (size_type i = 0; i < Q; ++i) {
+ gmm::sub_slice I2(i, ndof1/Q, Q);
+ gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
}
} else
- gmm::lu_inverse(M7);
+ { gmm::copy(M7, M7inv); gmm::lu_inverse(M7inv); }
- gmm::mult(M2, M1, MD);
+ 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(M7, M9, MPB);
+ 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(M7, MPC, MPD);
+ 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);
+ gmm::clean(M, 1E-13);
}
};
@@ -938,11 +948,11 @@ namespace getfem {
base_tensor t1, t2, ti, tv1, tv2, t1p, t2p;
base_matrix tv1p, tv2p, tvi;
- base_matrix M1(ndof2, ndof1), M2(ndof2, ndof2);
+ 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), M8(ndof1, ndof2);
+ base_matrix M7(ndof1, ndof1), M7inv(ndof1, ndof1), M8(ndof1, ndof2);
base_matrix M9(ndof1, ndof1), MD(ndof2, ndof1);
// Integrals on the element : \int_T Sym(Grad(D)).Grad(w) (M2)
@@ -997,12 +1007,12 @@ namespace getfem {
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof1; ++j)
- for (size_type k = 0; k < N*N; ++k)
+ for (size_type k = 0; k < N; ++k)
M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
- for (size_type k = 0; k < N*N; ++k)
+ for (size_type k = 0; k < N; ++k)
M8(i, j) += coeff * tv1p(i,k) * tv2p(j, k);
}
@@ -1051,68 +1061,68 @@ namespace getfem {
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof1; ++j)
- for (size_type k = 0; k < N*N; ++k)
+ for (size_type k = 0; k < N; ++k)
M7(i, j) += coeff * tv1p(i,k) * tv1p(j, k);
for (size_type i = 0; i < ndof1; ++i) // To be optimized
for (size_type j = 0; j < ndof2; ++j)
- for (size_type k = 0; k < N*N; ++k)
+ for (size_type k = 0; k < N; ++k)
M8(i, j) += coeff * tv1p(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 k = 0; k < N*N; ++k)
+ for (size_type k = 0; k < N; ++k)
M9(i, j) += coeff * tv1p(i,k) * tvi(j, k);
-
}
-
}
// Add the constraint with penalization
gmm::mult(gmm::transposed(M4), M4, aux2);
- gmm::add (aux2, M2);
+ gmm::add (gmm::scaled(aux2, 1E7), M2);
gmm::mult(gmm::transposed(M6), M6, aux2);
- gmm::add (aux2, M2);
+ gmm::add (gmm::scaled(aux2, 1E7), M2);
gmm::mult(gmm::transposed(M4), M3, aux1);
- gmm::add (aux1, M1);
+ gmm::add (gmm::scaled(aux1, 1E7), M1);
gmm::mult(gmm::transposed(M6), M5, aux1);
- gmm::add (aux1, M1);
+ gmm::add (gmm::scaled(aux1, 1E7), M1);
if (pf2->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof2, Q);
+ gmm::sub_slice I(0, ndof2/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- for (size_type i = 1; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof2, Q);
- gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2, I2, I2));
+ for (size_type i = 0; i < Q; ++i) {
+ gmm::sub_slice I2(i, ndof2/Q, Q);
+ gmm::copy(gmm::sub_matrix(M2, I, I), gmm::sub_matrix(M2inv, I2, I2));
}
} else
- gmm::lu_inverse(M2);
+ { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
if (pf1->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof1, Q);
+ gmm::sub_slice I(0, ndof1/Q, Q);
gmm::lu_inverse(gmm::sub_matrix(M7, I, I));
- for (size_type i = 1; i < Q; ++i) {
- gmm::sub_slice I2(i, ndof1, Q);
- gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7, I2, I2));
+ for (size_type i = 0; i < Q; ++i) {
+ gmm::sub_slice I2(i, ndof1/Q, Q);
+ gmm::copy(gmm::sub_matrix(M7, I, I), gmm::sub_matrix(M7inv, I2, I2));
}
} else
- gmm::lu_inverse(M7);
+ { gmm::copy(M7, M7inv); gmm::lu_inverse(M7inv); }
+
+ gmm::mult(M2inv, M1, MD);
+ gmm::clean(MD, gmm::vect_norminf(MD.as_vector()) * 1E-13);
- gmm::mult(M2, M1, MD);
-
// S = (I - inv(M7)*M9)(I - inv(M7)*M8*MD)
base_matrix MPB(ndof1, ndof1);
- gmm::mult(M7, M9, MPB);
+ 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(M7, MPC, MPD);
+ 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);
+ gmm::clean(M, 1E-13);
}
};
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index c432a66..38e1e04 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -4209,7 +4209,7 @@ namespace getfem {
void fem_precomp_pool::clear() {
for (const pfem_precomp &p : precomps)
- del_stored_object(p, true);
+ dal::del_stored_object(p, true);
precomps.clear();
}
diff --git a/src/getfem_fem_composite.cc b/src/getfem_fem_composite.cc
index 74d8187..74022f6 100644
--- a/src/getfem_fem_composite.cc
+++ b/src/getfem_fem_composite.cc
@@ -1,6 +1,6 @@
/*===========================================================================
- Copyright (C) 2002-2017 Yves Renard
+ Copyright (C) 2002-2019 Yves Renard
This file is a part of GetFEM++
@@ -122,6 +122,7 @@ namespace getfem {
std::vector<pdof_description> dofd(mf.nb_basic_dof());
for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
+ cout << "set poly of cv " << cv << endl;
pfem pf1 = mf.fem_of_element(cv);
if (!pf1->is_lagrange()) p->is_lagrange() = false;
if (!(pf1->is_equivalent())) p->is_equivalent() = false;
@@ -133,6 +134,7 @@ namespace getfem {
for (size_type k = 0; k < pf->nb_dof(cv); ++k) {
size_type igl = mf.ind_basic_dof_of_element(cv)[k];
base_poly fu = pf->base()[k];
+ cout << "adding " << fu << " : " << fu.dim() << endl;
base[igl].set_poly_of_subelt(cv, fu);
dofd[igl] = pf->dof_types()[k];
}
@@ -874,8 +876,8 @@ namespace getfem {
if (pf1 && (pf1->dim()+1 == pf0->dim()))
return phho->mf.fem_of_element(0);
}
-
- GMM_WARNING2("probably not a HHO method");
+
+ // GMM_WARNING2("probably not a HHO method");
return hho_method;
}
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 3614cd9..bf18df4 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -1142,7 +1142,7 @@ namespace getfem {
"transformation");
size_type ndof = Z.sizes()[0];
size_type Qmult = qdim / Z.sizes()[1];
- do_transformation(ndof*Qmult, t.sizes()[0]);
+ do_transformation(coeff_in.size(), ndof*Qmult);
return ga_instruction_val::exec();
}
@@ -1161,7 +1161,7 @@ namespace getfem {
GA_DEBUG_INFO("Instruction: gradient with elementary transformation");
size_type ndof = Z.sizes()[0];
size_type Qmult = qdim / Z.sizes()[1];
- do_transformation(ndof*Qmult, t.sizes()[0]);
+ do_transformation(coeff_in.size(), ndof*Qmult);
return ga_instruction_grad::exec();
}
@@ -1180,7 +1180,7 @@ namespace getfem {
GA_DEBUG_INFO("Instruction: Hessian with elementary transformation");
size_type ndof = Z.sizes()[0];
size_type Qmult = qdim / Z.sizes()[1];
- do_transformation(ndof*Qmult, t.sizes()[0]);
+ do_transformation(coeff_in.size(), ndof*Qmult);
return ga_instruction_hess::exec();
}
@@ -1199,7 +1199,7 @@ namespace getfem {
GA_DEBUG_INFO("Instruction: divergence with elementary transformation");
size_type ndof = Z.sizes()[0];
size_type Qmult = qdim / Z.sizes()[1];
- do_transformation(ndof*Qmult, t.sizes()[0]);
+ do_transformation(coeff_in.size(), ndof*Qmult);
return ga_instruction_diverg::exec();
}
@@ -1564,7 +1564,7 @@ namespace getfem {
size_type Qmult = qdim / Z.sizes()[1];
t_in.adjust_sizes(Qmult*ndof, Qmult*Z.sizes()[1]);
ga_instruction_copy_val_base::exec();
- do_transformation(ndof*Qmult, t_out.sizes()[0]);
+ do_transformation(t_out.sizes()[0], ndof*Qmult);
return 0;
}
@@ -1588,7 +1588,7 @@ namespace getfem {
size_type Qmult = qdim / Z.sizes()[1];
t_in.adjust_sizes(Qmult*ndof, Qmult*Z.sizes()[1], Z.sizes()[2]);
ga_instruction_copy_grad_base::exec();
- do_transformation(ndof*Qmult, t_out.sizes()[0]);
+ do_transformation(t_out.sizes()[0], ndof*Qmult);
return 0;
}
@@ -1612,7 +1612,7 @@ namespace getfem {
size_type Qmult = qdim / Z.sizes()[1];
t_in.adjust_sizes(Qmult*ndof, Qmult*Z.sizes()[1], Z.sizes()[2]);
ga_instruction_copy_hess_base::exec();
- do_transformation(ndof*Qmult, t_out.sizes()[0]);
+ do_transformation(t_out.sizes()[0], ndof*Qmult);
return 0;
}
@@ -1636,7 +1636,7 @@ namespace getfem {
size_type Qmult = qdim / Z.sizes()[1];
t_in.adjust_sizes(Qmult*ndof);
ga_instruction_copy_diverg_base::exec();
- do_transformation(ndof*Qmult, t_out.sizes()[0]);
+ do_transformation(t_out.sizes()[0], ndof*Qmult);
return 0;
}
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 719d4da..b211575 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -587,10 +587,10 @@ namespace getfem {
pnode->name = name;
if (ndt == 2) {
- std::string target_name = pnode->elementary_target;
- ga_parse_prefix_operator(target_name);
- ga_parse_prefix_test(target_name);
- pnode->elementary_target = target_name;
+ name = pnode->elementary_target;
+ ga_parse_prefix_operator(name);
+ ga_parse_prefix_test(name);
+ pnode->elementary_target = name;
}
// Group must be tested and it should be a fem variable
@@ -614,7 +614,7 @@ namespace getfem {
"Tensor with too many dimensions. Limited to 6");
if (test == 1) {
- pnode->name_test1 = name;
+ pnode->name_test1 = pnode->name;
pnode->interpolate_name_test1 = pnode->interpolate_name;
if (option == 1)
workspace.test1.insert
@@ -625,7 +625,7 @@ namespace getfem {
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
} else if (test == 2) {
- pnode->name_test2 = name;
+ pnode->name_test2 = pnode->name;
pnode->interpolate_name_test2 = pnode->interpolate_name;
if (option == 1)
workspace.test2.insert
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index e189e6e..13ea483 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -851,6 +851,23 @@ namespace getfem {
}
}
+ /* Select all the faces of the given mesh region (counted twice if they
+ are shared by two neighbour elements)
+ */
+ mesh_region all_faces_of_mesh(const mesh &m, const mesh_region &mr) {
+ mesh_region mrr;
+ mr.from_mesh(m);
+ mr.error_if_not_convexes();
+
+ for (mr_visitor i(mr); !i.finished(); ++i) {
+ size_type cv1 = i.cv();
+ short_type nbf = m.structure_of_convex(i.cv())->nb_faces();
+ for (short_type f = 0; f < nbf; ++f)
+ mrr.add(cv1, f);
+ }
+ return mrr;
+ }
+
/* Select all the faces sharing at least two element of the given mesh
region. Each face is represented only once and is arbitrarily chosen
between the two neighbour elements. Try to minimize the number of
- [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 <=
- [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