[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: |
Sun, 28 May 2017 11:22:52 -0400 (EDT) |
branch: devel-yves
commit 5a51ad5b1324eee76972234ec3f6b9aa82b2d174
Author: Yves Renard <address@hidden>
Date: Sun May 28 08:58:33 2017 +0200
A small optimization for interpolation
---
interface/src/Makefile.am | 1 -
interface/tests/matlab/demo_laplacian_pyramid.m | 114 ++++++++++++++++++++++++
interface/tests/python/Makefile.am | 1 +
src/bgeot_geometric_trans.cc | 56 ++++++++----
src/getfem/bgeot_poly.h | 5 ++
src/getfem/getfem_fem.h | 39 ++++++--
src/getfem_fem.cc | 12 ++-
7 files changed, 198 insertions(+), 30 deletions(-)
diff --git a/interface/src/Makefile.am b/interface/src/Makefile.am
index bf7dcce..9a858c5 100644
--- a/interface/src/Makefile.am
+++ b/interface/src/Makefile.am
@@ -25,7 +25,6 @@ endif
if BUILDSCILAB
subdirSCILAB=scilab
-.NOTPARALLEL: scilab
endif
SUBDIRS = . $(subdirMATLAB) $(subdirPYTHON) $(subdirSCILAB)
diff --git a/interface/tests/matlab/demo_laplacian_pyramid.m
b/interface/tests/matlab/demo_laplacian_pyramid.m
new file mode 100644
index 0000000..d3a6b9c
--- /dev/null
+++ b/interface/tests/matlab/demo_laplacian_pyramid.m
@@ -0,0 +1,114 @@
+% Copyright (C) 2005-2017 Yves Renard, Julien Pommier.
+%
+% This file is a part of GetFEM++
+%
+% GetFEM++ is free software; you can redistribute it and/or modify it
+% under the terms of the GNU Lesser General Public License as published
+% by the Free Software Foundation; either version 3 of the License, or
+% (at your option) any later version along with the GCC Runtime Library
+% Exception either version 3.1 or (at your option) any later version.
+% This program is distributed in the hope that it will be useful, but
+% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+% License and GCC Runtime Library Exception for more details.
+% You should have received a copy of the GNU Lesser General Public License
+% along with this program; if not, write to the Free Software Foundation,
+% Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+% Options for prescribing the Dirichlet condition
+dirichlet_version = 1; % 0 = simplification, 1 = with multipliers,
+ % 2 = penalization, 3 = Nitsche's method
+theta = 1; % Nitsche's method parameter theta
+gamma0 = 0.001; % Nitsche's method parameter gamma0 (gamma = gamma0*h)
+r = 1e8; % Penalization parameter
+draw = true;
+draw_mesh = false;
+
+K = 1; % Degree of the finite element method
+
+asize = size(who('automatic_var654'));
+if (asize(1)) draw = false; end;
+
+% trace on;
+gf_workspace('clear all');
+NX = 2;
+m = gf_mesh('pyramidal', [0:1/NX:1], [0:1/NX:1], [0:1/NX:1]);
+
+
+% Create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
+mf = gf_mesh_fem(m,1);
+% Assign the pyramidal fem to all convexes of the mesh_fem, and define an
+% integration method
+gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_PYRAMID_LAGRANGE(%d)', K)));
+mim = gf_mesh_im(m, gf_integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(5))'));
+
+
+% Detect the border of the mesh
+border = gf_mesh_get(m, 'outer faces');
+% Mark it as boundary GAMMAD=1
+GAMMAD=1;
+gf_mesh_set(m, 'region', GAMMAD, border);
+if (draw_mesh)
+ figure(1);
+ gf_plot_mesh(m, 'regions', [GAMMAD], 'convexes', 'on'); % the boundary edges
appear in red
+ pause(1);
+end
+
+% Interpolate the exact solution
+% Uexact = gf_mesh_fem_get(mf, 'eval', { '10*y.*(y-1).*x.*(x-1)+10*x.^5' });
+Uexact = gf_mesh_fem_get(mf, 'eval', { 'x.*sin(2*pi*x).*sin(2*pi*y)' });
+% Opposite of its laplacian
+% F = gf_mesh_fem_get(mf, 'eval', {
'-(20*(x.^2+y.^2)-20*x-20*y+200*x.^3)' });
+F = gf_mesh_fem_get(mf, 'eval', { '4*pi*(2*pi*x.*sin(2*pi*x) -
cos(2*pi*x)).*sin(2*pi*y)' });
+
+md=gf_model('real');
+gf_model_set(md, 'add fem variable', 'u', mf);
+gf_model_set(md, 'add linear generic assembly brick', mim,
'Grad_u.Grad_Test_u');
+% gf_model_set(md, 'add Laplacian brick', mim, 'u');
+gf_model_set(md, 'add initialized fem data', 'VolumicData', mf, F);
+gf_model_set(md, 'add source term brick', mim, 'u', 'VolumicData');
+gf_model_set(md, 'add initialized fem data', 'DirichletData', mf, Uexact);
+switch (dirichlet_version)
+ case 0,
+ gf_model_set(md, 'add Dirichlet condition with simplification', 'u',
GAMMAD, 'DirichletData');
+ case 1,
+ gf_model_set(md, 'add Dirichlet condition with multipliers', mim, 'u', mf,
GAMMAD, 'DirichletData');
+ case 2,
+ gf_model_set(md, 'add Dirichlet condition with penalization', mim, 'u', r,
GAMMAD, 'DirichletData');
+ case 3,
+ gf_model_set(md, 'add initialized data', 'gamma0', [gamma0]);
+ expr = gf_model_get(md, 'Neumann term', 'u', GAMMAD);
+ gf_model_set(md, 'add Dirichlet condition with Nitsche method', mim, 'u',
expr, 'gamma0', GAMMAD, theta, 'DirichletData');
+end
+gf_model_get(md, 'solve');
+U = gf_model_get(md, 'variable', 'u');
+Utest = U*0;
+if (draw)
+ figure(2);
+ subplot(2,1,1); gf_plot(mf,U,'mesh','off', 'cvlst', gf_mesh_get(m, 'outer
faces'), 'refine', 4);
+ colorbar; title('computed solution');
+
+ subplot(2,1,2); gf_plot(mf,Uexact,'mesh','on', 'cvlst', gf_mesh_get(m,
'outer faces'), 'refine', 4);
+ colorbar;title('exact solution');
+end
+
+if(0)
+ for i = 1:size(U,2)
+ Utest(i) = 1;
+ figure(3);
+ gf_plot(mf,Utest,'mesh','on', 'cvlst', gf_mesh_get(m, 'outer faces'),
'refine', 8);
+ colorbar;title('shape function');
+ Utest(i) = 0;
+ pause;
+ end
+end
+
+err = gf_compute(mf, Uexact-U, 'H1 norm', mim);
+
+disp(sprintf('H1 norm of error: %g', err));
+
+if (err > 2E-4)
+ error('Laplacian test: error to big');
+end
+
+
diff --git a/interface/tests/python/Makefile.am
b/interface/tests/python/Makefile.am
index 6caa3c8..f80d0a9 100644
--- a/interface/tests/python/Makefile.am
+++ b/interface/tests/python/Makefile.am
@@ -22,6 +22,7 @@ EXTRA_DIST= \
demo_crack.py \
demo_fictitious_domains.py \
demo_laplacian.py \
+ demo_laplacian_pyramid.py \
demo_laplacian_DG.py \
demo_laplacian_aposteriori.py \
demo_mortar.py \
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 7f1379d..8ff4835 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -483,14 +483,43 @@ namespace bgeot {
struct igeometric_trans : public geometric_trans {
std::vector<FUNC> trans;
+ mutable std::vector<std::vector<FUNC>> grad_, hess_;
+
+ void compute_grad_() const {
+ size_type R = trans.size();
+ dim_type n = dim();
+ grad_.resize(R);
+ for (size_type i = 0; i < R; ++i) {
+ grad_[i].resize(n);
+ for (dim_type j = 0; j < n; ++j) {
+ grad_[i][j] = trans[i]; grad_[i][j].derivative(j);
+ }
+ }
+ }
+ void compute_hess_() const {
+ size_type R = trans.size();
+ dim_type n = dim();
+ hess_.resize(R);
+ for (size_type i = 0; i < R; ++i) {
+ hess_[i].resize(n*n);
+ for (dim_type j = 0; j < n; ++j) {
+ for (dim_type k = 0; k < n; ++k) {
+ hess_[i][j+k*n] = trans[i];
+ hess_[i][j+k*n].derivative(j); hess_[i][j+k*n].derivative(k);
+ }
+ }
+ }
+ }
+
virtual void poly_vector_val(const base_node &pt, base_vector &val) const {
val.resize(nb_points());
for (size_type k = 0; k < nb_points(); ++k)
val[k] = to_scalar(trans[k].eval(pt.begin()));
}
- virtual void poly_vector_val(const base_node &pt, const convex_ind_ct
&ind_ct,
+ virtual void poly_vector_val(const base_node &pt,
+ const convex_ind_ct &ind_ct,
base_vector &val) const {
size_type nb_funcs=ind_ct.size();
val.resize(nb_funcs);
@@ -499,40 +528,35 @@ namespace bgeot {
}
virtual void poly_vector_grad(const base_node &pt, base_matrix &pc) const {
+ if (!(grad_.size())) compute_grad_();
FUNC PP;
pc.base_resize(nb_points(),dim());
for (size_type i = 0; i < nb_points(); ++i)
- for (dim_type n = 0; n < dim(); ++n) {
- PP = trans[i];
- PP.derivative(n);
- pc(i, n) = to_scalar(PP.eval(pt.begin()));
- }
+ for (dim_type n = 0; n < dim(); ++n)
+ pc(i, n) = to_scalar(grad_[i][n].eval(pt.begin()));
}
virtual void poly_vector_grad(const base_node &pt,
const convex_ind_ct &ind_ct,
base_matrix &pc) const {
+ if (!(grad_.size())) compute_grad_();
FUNC PP;
size_type nb_funcs=ind_ct.size();
pc.base_resize(nb_funcs,dim());
for (size_type i = 0; i < nb_funcs; ++i)
- for (dim_type n = 0; n < dim(); ++n) {
- PP = trans[ind_ct[i]];
- PP.derivative(n);
- pc(i, n) = to_scalar(PP.eval(pt.begin()));
- }
+ for (dim_type n = 0; n < dim(); ++n)
+ pc(i, n) = to_scalar(grad_[ind_ct[i]][n].eval(pt.begin()));
}
virtual void poly_vector_hess(const base_node &pt, base_matrix &pc) const {
+ if (!(grad_.size())) compute_grad_();
FUNC PP, QP;
pc.base_resize(nb_points(),dim()*dim());
for (size_type i = 0; i < nb_points(); ++i)
for (dim_type n = 0; n < dim(); ++n) {
- QP = trans[i]; QP.derivative(n);
- for (dim_type m = 0; m <= n; ++m) {
- PP = QP; PP.derivative(m);
- pc(i, n*dim()+m) = pc(i, m*dim()+n) =
to_scalar(PP.eval(pt.begin()));
- }
+ for (dim_type m = 0; m <= n; ++m)
+ pc(i, n*dim()+m) = pc(i, m*dim()+n) =
+ to_scalar(hess_[i][m*dim()+n].eval(pt.begin()));
}
}
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index 1f8112f..3c4275c 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -749,6 +749,11 @@ namespace bgeot
return R;
}
+ template<typename T> std::ostream &operator <<
+ (std::ostream &o, const rational_fraction<T>& P) {
+ o << "[" << P.numerator() << "]/[" << P.denominator() << "]";
+ return o;
+ }
typedef rational_fraction<opt_long_scalar_type> base_rational_fraction;
diff --git a/src/getfem/getfem_fem.h b/src/getfem/getfem_fem.h
index 115b10f..56721a7 100644
--- a/src/getfem/getfem_fem.h
+++ b/src/getfem/getfem_fem.h
@@ -485,7 +485,35 @@ namespace getfem {
template <class FUNC> class fem : public virtual_fem {
protected :
std::vector<FUNC> base_;
+ mutable std::vector<std::vector<FUNC>> grad_, hess_;
+ void compute_grad_() const {
+ size_type R = nb_base_components(0);
+ dim_type n = dim();
+ grad_.resize(R);
+ for (size_type i = 0; i < R; ++i) {
+ grad_[i].resize(n);
+ for (dim_type j = 0; j < n; ++j) {
+ grad_[i][j] = base_[i]; grad_[i][j].derivative(j);
+ }
+ }
+ }
+
+ void compute_hess_() const {
+ size_type R = nb_base_components(0);
+ dim_type n = dim();
+ hess_.resize(R);
+ for (size_type i = 0; i < R; ++i) {
+ hess_[i].resize(n*n);
+ for (dim_type j = 0; j < n; ++j) {
+ for (dim_type k = 0; k < n; ++k) {
+ hess_[i][j+k*n] = base_[i];
+ hess_[i][j+k*n].derivative(j); hess_[i][j+k*n].derivative(k);
+ }
+ }
+ }
+ }
+
public :
/// Gives the array of basic functions (components).
@@ -506,6 +534,7 @@ namespace getfem {
reference element directions 0,..,dim-1 and returns the result in
t(nb_base,target_dim,dim) */
void grad_base_value(const base_node &x, base_tensor &t) const {
+ if (!(grad_.size())) compute_grad_();
bgeot::multi_index mi(3);
dim_type n = dim();
mi[2] = n; mi[1] = target_dim(); mi[0] = short_type(nb_base(0));
@@ -514,12 +543,13 @@ namespace getfem {
base_tensor::iterator it = t.begin();
for (dim_type j = 0; j < n; ++j)
for (size_type i = 0; i < R; ++i, ++it)
- { FUNC f = base_[i]; f.derivative(j); *it =
bgeot::to_scalar(f.eval(x.begin())); }
+ *it = bgeot::to_scalar(grad_[i][j].eval(x.begin()));
}
/** Evaluates at point x, the hessian of all base functions w.r.t. the
reference element directions 0,..,dim-1 and returns the result in
t(nb_base,target_dim,dim,dim) */
void hess_base_value(const base_node &x, base_tensor &t) const {
+ if (!(hess_.size())) compute_hess_();
bgeot::multi_index mi(4);
dim_type n = dim();
mi[3] = n; mi[2] = n; mi[1] = target_dim();
@@ -529,11 +559,8 @@ namespace getfem {
base_tensor::iterator it = t.begin();
for (dim_type k = 0; k < n; ++k)
for (dim_type j = 0; j < n; ++j)
- for (size_type i = 0; i < R; ++i, ++it) {
- FUNC f = base_[i];
- f.derivative(j); f.derivative(k);
- *it = bgeot::to_scalar(f.eval(x.begin()));
- }
+ for (size_type i = 0; i < R; ++i, ++it)
+ *it = bgeot::to_scalar(hess_[i][j+k*n].eval(x.begin()));
}
};
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index efcc1c2..888cf7b 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1280,11 +1280,11 @@ namespace getfem {
} else if (k == 2) {
p->base().resize(14);
- base_poly xi0 = bgeot::read_base_poly(3, "(1-z-x)*0.5");
- base_poly xi1 = bgeot::read_base_poly(3, "(1-z-y)*0.5");
- base_poly xi2 = bgeot::read_base_poly(3, "(1-z+x)*0.5");
- base_poly xi3 = bgeot::read_base_poly(3, "(1-z+y)*0.5");
- base_poly z = bgeot::read_base_poly(3, "z");
+ base_poly xi0 = bgeot::read_base_poly(3, "(1-z-x)*0.5");
+ base_poly xi1 = bgeot::read_base_poly(3, "(1-z-y)*0.5");
+ base_poly xi2 = bgeot::read_base_poly(3, "(1-z+x)*0.5");
+ base_poly xi3 = bgeot::read_base_poly(3, "(1-z+y)*0.5");
+ base_poly z = bgeot::read_base_poly(3, "z");
base_poly un_z = bgeot::read_base_poly(3, "1-z");
bgeot::base_rational_fraction Q(bgeot::read_base_poly(3, "1"), un_z);
@@ -1318,8 +1318,6 @@ namespace getfem {
p->add_node(lag_dof, base_small_vector( 0.5, 0.5, 0.5));
p->add_node(lag_dof, base_small_vector( 0.0, 0.0, 1.0));
- GMM_ASSERT1(false, "to be completed");
-
} else GMM_ASSERT1(false, "Sorry, pyramidal Lagrange fem "
"implemented only for degree 0, 1 or 2");