[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 43078bff190a072beb91acce5a32f4c95f457879
Author: Yves Renard <address@hidden>
Date: Wed Aug 21 16:41:29 2019 +0200
add a test on elasticity problem for HHO method
---
.gitignore | 1 +
doc/sphinx/source/userdoc/hho.rst | 12 +-
interface/tests/python/Makefile.am | 3 +
interface/tests/python/demo_elasticity_HHO.py | 182 ++++++++++++++++++++++++++
interface/tests/python/demo_laplacian_HHO.py | 48 +++----
interface/tests/python/demo_wave_equation.py | 4 +-
src/getfem_HHO.cc | 46 +++----
7 files changed, 224 insertions(+), 72 deletions(-)
diff --git a/.gitignore b/.gitignore
index b7a5bfb..3123daf 100644
--- a/.gitignore
+++ b/.gitignore
@@ -129,6 +129,7 @@ interface/tests/python/results1/
/doc/sphinx/source/python/cmdref*
/doc/sphinx/source/matlab/cmdref*
/doc/sphinx/build/
+/doc/sphinx/homepage/
/doc/sphinx/*.pdf
/doc/sphinx/source/matlab/images/hierarchy.png
/doc/sphinx/source/project/images/diagram.png
diff --git a/doc/sphinx/source/userdoc/hho.rst
b/doc/sphinx/source/userdoc/hho.rst
index b7ce053..a5a3086 100644
--- a/doc/sphinx/source/userdoc/hho.rst
+++ b/doc/sphinx/source/userdoc/hho.rst
@@ -45,9 +45,9 @@ In order to be used, the elementary transformation
corresponding to this operato
add_HHO_reconstructed_gradient(model, transname);
-where ``transname`` is an arbitrary name which will designate the
transformation in the weak form language. Then, it will be possible to refer to
the reconstructed gradient of a variable ``u`` into the weak form language as
``Elementary_transformation(u, HHO_grad, ur)``, if ``transname="HHO_grad"``.
The third parameter of the transformation ``ur`` should be a fem variable or a
data of the model. This variable will not be used on itself but will determine
the finite element space of the r [...]
+where ``transname`` is an arbitrary name which will designate the
transformation in the weak form language. Then, it will be possible to refer to
the reconstructed gradient of a variable ``u`` into the weak form language as
``Elementary_transformation(u, HHO_grad, Gu)``, if ``transname="HHO_grad"``.
The third parameter of the transformation ``Gu`` should be a fem variable or a
data of the model. This variable will not be used on itself but will determine
the finite element space of the r [...]
-This is an example of use in Python for a two-dimensional mesh ``m`` ::
+This is an example of use with the Python interface for a two-dimensional
triangule mesh ``m`` ::
mfu = gf.MeshFem(m, 1)
mfgu = gf.MeshFem(m, N)
@@ -66,7 +66,7 @@ The macro definitions allowing to use the gradient of the
variable inside weak f
md.add_linear_term(mim, 'HHO_Grad_u.HHO_Grad_Test_u')
-A complete example of use is given in the test program
:file:`interface/tests/demo_laplacian_HHO.py`
+Two complete examples of use are given in the test programs
:file:`interface/tests/demo_laplacian_HHO.py` and
:file:`interface/tests/demo_elasticity_HHO.py`.
Reconstructed symmetrized gradient
++++++++++++++++++++++++++++++++++
@@ -82,7 +82,7 @@ The elementary transformation corresponding to this operator
can be added to the
add_HHO_reconstructed_symmetrized_gradient(model, transname);
-and then be used into the weak form language as ``Elementary_transformation(u,
HHO_sym_grad, ur)``, if ``transname="HHO_sym_grad"``, with ``ur`` still
determining the reconstruction space.
+and then be used into the weak form language as ``Elementary_transformation(u,
HHO_sym_grad, Gu)``, if ``transname="HHO_sym_grad"``, with ``Gu`` still
determining the reconstruction space.
Reconstructed variable
++++++++++++++++++++++
@@ -144,9 +144,9 @@ For vector field variables having the same dimension as the
domain, there exists
.. math::
S^s(u) = \Pi_{\partial T}(u_{\partial T} - D^s(u) - \Pi_{T}(u_T - D^s(u)))
-The corresponding elementary transformation can be added to the model by the
two commands::
+The corresponding elementary transformations can be added to the model by the
two commands::
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. A example
of use is also given in the test program
:file:`interface/tests/demo_laplacian_HHO.py`.
+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`.
diff --git a/interface/tests/python/Makefile.am
b/interface/tests/python/Makefile.am
index 113a191..e0f42d8 100644
--- a/interface/tests/python/Makefile.am
+++ b/interface/tests/python/Makefile.am
@@ -52,6 +52,8 @@ EXTRA_DIST= \
demo_thermo_elasticity_electrical_coupling.py \
demo_cracked_thermo_elastic_body.py \
demo_nonlinear_elasticity.py \
+ demo_elasticity_HHO.py \
+ demo_elastic_ring_contact.py \
demo_wheel_contact.py \
demo_tripod.py \
demo_tripod_alt.py \
@@ -70,6 +72,7 @@ TESTS = \
demo_wave_equation.py \
demo_laplacian.py \
demo_laplacian_HHO.py \
+ demo_elasticity_HHO.py \
$(optpy)
AM_TESTS_ENVIRONMENT = \
diff --git a/interface/tests/python/demo_elasticity_HHO.py
b/interface/tests/python/demo_elasticity_HHO.py
new file mode 100644
index 0000000..c3653d9
--- /dev/null
+++ b/interface/tests/python/demo_elasticity_HHO.py
@@ -0,0 +1,182 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2019-2019 Yves Renard.
+#
+# 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 2.1 of the License, 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 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.
+#
+############################################################################
+""" 2D elasticity 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++.
+
+ $Id$
+"""
+# Import basic modules
+import getfem as gf
+import numpy as np
+
+## Parameters
+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_symmetric_gradient = True # Use symmetric gradient reconstruction or
not
+
+E = 1 # Young's modulus
+nu = 0.499 # Poisson ratio
+
+cmu = E/(2*(1+nu)) # Lame coefficient
+clambda = 2*cmu*nu/(1-2*nu) # Lame coefficient
+
+# 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];')
+
+N = m.dim();
+
+# Meshfems
+mfu = gf.MeshFem(m, N)
+mfgu = gf.MeshFem(m, N, N)
+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)'))
+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)'))
+
+print('nbdof : %d' % mfu.nbdof());
+
+# Integration method used
+mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
+
+# Boundary selection
+flst = m.outer_faces()
+GAMMAD = 1
+m.set_region(GAMMAD, flst)
+
+# 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)
+a = 8.
+Ue = mfur.eval('[np.sin((%g)*x), -(%g)*y*np.cos((%g)*x)]' % (a,a,a),
globals(), locals())
+
+# Interpolate the source term
+F1 = mfrhs.eval('(%g)*(pow(%g,2.))*np.sin((%g)*x),
-(%g)*(pow(%g,3.))*y*np.cos((%g)*x)' % (cmu, a, a, cmu, a, a), globals(),
locals())
+
+
+
+# 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):
+ if (using_symmetric_gradient):
+ md.add_HHO_reconstructed_symmetrized_gradient('HHO_Grad')
+ md.add_HHO_reconstructed_gradient('HHO_vGrad')
+ md.add_HHO_reconstructed_symmetrized_value('HHO_Val')
+ md.add_HHO_symmetrized_stabilization('HHO_Stab')
+ md.add_macro('HHO_vGrad_u', 'Elementary_transformation(u, HHO_vGrad, Gu)')
+ else :
+ md.add_HHO_reconstructed_gradient('HHO_Grad')
+ md.add_HHO_reconstructed_value('HHO_Val')
+ md.add_HHO_stabilization('HHO_Stab')
+ md.add_macro('HHO_vGrad_u', 'Elementary_transformation(u, HHO_Grad, Gu)')
+
+
+ md.add_macro('HHO_Val_u', 'Elementary_transformation(u, HHO_Val, ur)')
+ 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)')
+
+
+# Elasticity term on u
+md.add_initialized_data('cmu', [cmu])
+md.add_initialized_data('clambda', [clambda])
+if (using_HHO):
+ # Elasticity term
+ md.add_linear_term(mim, 'clambda*Trace(HHO_Grad_u)*Trace(HHO_Grad_Test_u)'
+ + '+ 2*cmu*Sym(HHO_Grad_u):Sym(HHO_Grad_Test_u)')
+ # Stabilization term
+ md.add_linear_term(mim, 'cmu*HHO_Stab_u.HHO_Stab_Test_u', ALL_FACES)
+else:
+ md.add_isotropic_linearized_elasticity_brick(mim, 'u', 'clambda', 'cmu')
+
+# Volumic source term
+md.add_initialized_fem_data('VolumicData', mfrhs, F1)
+md.add_source_term_brick(mim, 'u', 'VolumicData')
+
+# Dirichlet condition
+md.add_initialized_fem_data("Ue", mfur, Ue)
+
+if (Dirichlet_with_multipliers):
+ md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu, GAMMAD, 'Ue')
+else:
+ md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
+ GAMMAD, 'Ue')
+
+# Assembly of the linear system and solve.
+md.solve()
+
+# Error computation
+U = md.variable('u')
+L2error = gf.asm('generic', mim, 0, 'Norm_sqr(u-Ue)', -1, md)
+H1error = gf.asm('generic', mim, 0, 'Norm_sqr(Grad_u-Grad_Ue)', -1, md)
+H1error = np.sqrt(L2error + H1error); L2error = np.sqrt(L2error)
+print('Error in L2 norm (without reconstruction): %g' % L2error)
+print('Error in H1 norm (without reconstruction): %g' % H1error)
+
+if (using_HHO):
+ L2error = gf.asm('generic', mim, 0, 'Norm_sqr(HHO_Val_u-Ue)', -1, md)
+ H1error = gf.asm('generic', mim, 0, 'Norm_sqr(HHO_vGrad_u-Grad_Ue)', -1, md)
+ H1error = np.sqrt(L2error + H1error); L2error = np.sqrt(L2error)
+ print('Error in L2 norm (with reconstruction): %g' % L2error)
+ print('Error in H1 norm (with reconstruction): %g' % H1error)
+
+
+# Export data
+# mfur.export_to_pos('elasticity_e.pos', Ue, 'Exact solution')
+mfu.export_to_pos('elasticity.pos', U, 'Computed solution')
+print('You can view the solution with (for example):')
+print('gmsh elasticity.pos')
+
+
+if (H1error > 0.013):
+ print('Error too large !')
+ exit(1)
diff --git a/interface/tests/python/demo_laplacian_HHO.py
b/interface/tests/python/demo_laplacian_HHO.py
index 06b472d..db8de14 100644
--- a/interface/tests/python/demo_laplacian_HHO.py
+++ b/interface/tests/python/demo_laplacian_HHO.py
@@ -31,7 +31,7 @@ import getfem as gf
import numpy as np
## Parameters
-NX = 40 # Mesh parameter.
+NX = 10 # Mesh parameter.
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
@@ -135,59 +135,41 @@ md.add_normal_source_term_brick(mim, 'u', 'NeumannData',
NEUMANN_BOUNDARY_NUM)
# Dirichlet condition on the left.
-md.add_initialized_fem_data("DirichletData", mfur, Ue)
+md.add_initialized_fem_data("Ue", mfur, Ue)
if (Dirichlet_with_multipliers):
md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu,
- DIRICHLET_BOUNDARY_NUM1,
- 'DirichletData')
+ DIRICHLET_BOUNDARY_NUM1, 'Ue')
else:
md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
- DIRICHLET_BOUNDARY_NUM1,
- 'DirichletData')
+ DIRICHLET_BOUNDARY_NUM1, 'Ue')
# Dirichlet condition on the top.
# Two Dirichlet brick in order to test the multiplier
# selection in the intersection.
if (Dirichlet_with_multipliers):
md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu,
- DIRICHLET_BOUNDARY_NUM2,
- 'DirichletData')
+ DIRICHLET_BOUNDARY_NUM2, 'Ue')
else:
md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
- DIRICHLET_BOUNDARY_NUM2,
- 'DirichletData')
-# gf.memstats()
-# md.listvar()
-# md.listbricks()
+ DIRICHLET_BOUNDARY_NUM2, 'Ue')
# Assembly of the linear system and solve.
md.solve()
-# Main unknown
+# Error computation
U = md.variable('u')
+L2error = gf.asm('generic', mim, 0, 'sqr(u-Ue)', -1, md)
+H1error = gf.asm('generic', mim, 0, 'Norm_sqr(Grad_u-Grad_Ue)', -1, md)
+H1error = np.sqrt(L2error + H1error); L2error = np.sqrt(L2error)
+print('Error in L2 norm (without reconstruction): %g' % L2error)
+print('Error in H1 norm (without reconstruction): %g' % H1error)
if (using_HHO):
- L2error = gf.asm('generic', mim, 0, 'sqr(u-Ue)',
- -1, md, 'Ue', 0, mfur, Ue)
- H1error = gf.asm('generic', mim, 0, 'Norm_sqr(Grad_u-Grad_Ue)',
- -1, md, 'Ue', 0, mfur, Ue)
- H1error = np.sqrt(L2error + H1error)
- L2error = np.sqrt(L2error)
- print('Error in L2 norm (without reconstruction): %g' % L2error)
- print('Error in H1 norm (without reconstruction): %g' % H1error)
- L2error = gf.asm('generic', mim, 0, 'sqr(HHO_Val_u-Ue)',
- -1, md, 'Ue', 0, mfur, Ue)
- H1error = gf.asm('generic', mim, 0, 'Norm_sqr(HHO_Grad_u-Grad_Ue)',
- -1, md, 'Ue', 0, mfur, Ue)
- H1error = np.sqrt(L2error + H1error)
- L2error = np.sqrt(L2error)
+ L2error = gf.asm('generic', mim, 0, 'sqr(HHO_Val_u-Ue)', -1, md)
+ H1error = gf.asm('generic', mim, 0, 'Norm_sqr(HHO_Grad_u-Grad_Ue)', -1, md)
+ H1error = np.sqrt(L2error + H1error); L2error = np.sqrt(L2error)
print('Error in L2 norm (with reconstruction): %g' % L2error)
print('Error in H1 norm (with reconstruction): %g' % H1error)
-else :
- L2error = gf.compute(mfu, U-Ue, 'L2 norm', mim)
- H1error = gf.compute(mfu, U-Ue, 'H1 norm', mim)
- print('Error in L2 norm : %g' % L2error)
- print('Error in H1 norm : %g' % H1error)
# Export data
diff --git a/interface/tests/python/demo_wave_equation.py
b/interface/tests/python/demo_wave_equation.py
index 1deb5a1..693f1d3 100644
--- a/interface/tests/python/demo_wave_equation.py
+++ b/interface/tests/python/demo_wave_equation.py
@@ -1,5 +1,5 @@
#!/usr/bin/env python
-# -*- coding: UTF8 -*-
+# -*- coding: utf-8 -*-
# Python GetFEM++ interface
#
# Copyright (C) 2015-2019 FABRE Mathieu, SECK Mamadou, DALLERIT Valentin,
@@ -31,7 +31,7 @@ import numpy as np
import getfem as gf
import os
-NX = 20
+NX = 10
m = gf.Mesh('cartesian', np.arange(0., 1.+1./NX,1./NX),
np.arange(0., 1.+1./NX,1./NX));
diff --git a/src/getfem_HHO.cc b/src/getfem_HHO.cc
index c9d26a8..7831047 100644
--- a/src/getfem_HHO.cc
+++ b/src/getfem_HHO.cc
@@ -428,7 +428,7 @@ namespace getfem {
scalar_type b(0), a = coeff *
(tv1p(i, k1) - (i < ndofi ? tvi(i, k1) : 0.));
for (size_type k2 = 0; k2 < N; ++k2)
- b += a * tv2.as_vector()[j+(k1 + k2*Q)*ndof2] * un[k2];
+ b += a * tv2.as_vector()[j+(k1+k2*Q)*ndof2] * un[k2];
M1(j, i) += b;
}
}
@@ -578,7 +578,7 @@ namespace getfem {
}
// Integrals on the faces : \int_{dT}(v_{dT} - v_T).(Sym(Grad(w)).n) (M1)
- // \int_{dT} Skew(n x v_{dT} - v_{dT} x n) (M5)
+ // \int_{dT} 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);
size_type first_ind = pim->ind_first_point_on_face(ifc);
@@ -605,16 +605,16 @@ namespace getfem {
scalar_type b(0), a = coeff *
(tv1p(i, k1) - (i < ndofi ? tvi(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];
+ b += a * 0.5 * (tv2.as_vector()[j+(k1+k2*N)*ndof2] +
+ tv2.as_vector()[j+(k2+k1*N)*ndof2]) * un[k2];
M1(j, i) += b;
}
for (size_type i = 0; i < ndof1; ++i)
for (size_type k1 = 0; k1 < N; ++k1)
for (size_type k2 = 0; k2 < N; ++k2)
- M5(k1+k2*N, i) += 0.5 * coeff * (tv1p(i, k2) * un[k1] -
- tv1p(i, k1) * un[k2]);
+ M5(k1+k2*N, i) += 0.5 * coeff * (tv1p(i, k1) * un[k2] -
+ tv1p(i, k2) * un[k1]);
}
}
@@ -622,21 +622,13 @@ namespace getfem {
gmm::mult(gmm::transposed(M4), M4, aux2);
gmm::add (gmm::scaled(aux2, 1.E7), M2);
gmm::mult(gmm::transposed(M6), M6, aux2);
- gmm::add (gmm::scaled(aux2, 1.E7), M2);
+ gmm::add (gmm::scaled(aux2, 1.E5), M2);
gmm::mult(gmm::transposed(M4), M3, aux1);
gmm::add (gmm::scaled(aux1, 1.E7), M1);
gmm::mult(gmm::transposed(M6), M5, aux1);
- gmm::add (gmm::scaled(aux1, 1.E7), M1);
+ gmm::add (gmm::scaled(aux1, 1.E5), M1);
- if (pf2->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof2/Q, Q);
- gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- 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::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
+ gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv);
gmm::mult(M2inv, M1, M);
gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
@@ -971,7 +963,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)
@@ -1054,8 +1046,8 @@ namespace getfem {
for (size_type i = 0; i < ndof1; ++i)
for (size_type k1 = 0; k1 < N; ++k1)
for (size_type k2 = 0; k2 < N; ++k2)
- M5(k1+k2*N, i) += 0.5 * coeff * (tv1p(i, k2) * un[k1] -
- tv1p(i, k1) * un[k2]);
+ 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)
@@ -1078,21 +1070,13 @@ namespace getfem {
gmm::mult(gmm::transposed(M4), M4, aux2);
gmm::add (gmm::scaled(aux2, 1E7), M2);
gmm::mult(gmm::transposed(M6), M6, aux2);
- gmm::add (gmm::scaled(aux2, 1E7), M2);
+ gmm::add (gmm::scaled(aux2, 1E5), M2);
gmm::mult(gmm::transposed(M4), M3, aux1);
gmm::add (gmm::scaled(aux1, 1E7), M1);
gmm::mult(gmm::transposed(M6), M5, aux1);
- gmm::add (gmm::scaled(aux1, 1E7), M1);
+ gmm::add (gmm::scaled(aux1, 1E5), M1);
- if (pf2->target_dim() == 1 && Q > 1) {
- gmm::sub_slice I(0, ndof2/Q, Q);
- gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
- 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::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
+ gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv);
if (pf1->target_dim() == 1 && Q > 1) {
gmm::sub_slice I(0, ndof1/Q, Q);
- [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 <=
- [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, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27