[getfem] 02/05: Merge.
Anton Gladky
gladk at alioth.debian.org
Sat Nov 2 20:16:59 UTC 2013
This is an automated email from the git hooks/post-receive script.
gladk pushed a commit to branch master
in repository getfem.
commit 5fabc59d98361da1b9d9918fc6c35fddd830bc38
Merge: fa2b872 e3bfac6
Author: Anton Gladky <gladky.anton at gmail.com>
Date: Sat Nov 2 19:44:14 2013 +0100
Merge.
bin/extract_doc | 11 +-
configure.ac | 52 +-
contrib/level_set_contact/Makefile.am | 4 +-
contrib/level_set_contact/contact_problem.h | 4 +-
contrib/xfem_large_strain/vertex_large_strain.cc | 9 +
doc/sphinx/source/.templates/indexcontent.html | 2 +-
doc/sphinx/source/biblio.rst | 35 +-
doc/sphinx/source/python/install.rst | 2 +
doc/sphinx/source/python/pygf.rst | 6 +
.../source/userdoc/images/getfemusercorrection.png | Bin 10640 -> 7608 bytes
.../source/userdoc/images/getfemusertransition.png | Bin 0 -> 11930 bytes
doc/sphinx/source/userdoc/model_continuation.rst | 211 ++--
doc/sphinx/source/userdoc/model_solvers.rst | 28 +-
doc/sphinx/source/userdoc/parallel.rst | 130 ++-
interface/src/getfemint.h | 2 +-
interface/src/gf_compute.cc | 34 +-
interface/src/gf_cont_struct.cc | 5 +
interface/src/gf_cont_struct_get.cc | 7 +-
interface/src/gf_mesh_fem_get.cc | 27 +-
interface/src/gf_model_get.cc | 2 +-
interface/src/matlab/Makefile.am | 2 +-
interface/src/python/Makefile.am | 17 +-
interface/src/scilab/Makefile.am | 3 +
interface/src/scilab/demos/data/vee_h_0.03.mesh | 906 +++++++++++++++
interface/src/scilab/demos/demo_continuation.sce | 11 +-
.../src/scilab/demos/demo_continuation_block.sce | 349 ++++++
.../src/scilab/demos/demo_continuation_vee.sce | 352 ++++++
interface/tests/matlab/demo_dynamic_contact.m | 17 +-
interface/tests/matlab/demo_static_contact.m | 14 +-
interface/tests/python/demo_laplacian.py | 4 +-
interface/tests/python/demo_parallel_laplacian.py | 54 +-
src/Makefile.am | 1 +
src/dal_static_stored_objects.cc | 210 ++--
src/getfem/bgeot_tensor.h | 246 +++--
src/getfem/getfem_assembling_tensors.h | 58 +-
src/getfem/getfem_config.h | 123 ++-
src/getfem/getfem_continuation.h | 26 +-
src/getfem/getfem_convect.h | 28 +-
src/getfem/getfem_derivatives.h | 14 +-
src/getfem/getfem_error_estimate.h | 18 +-
src/getfem/getfem_im_list.h | 17 +-
src/getfem/getfem_interpolation.h | 7 +-
src/getfem/getfem_mesh_region.h | 16 +-
src/getfem/getfem_model_solvers.h | 218 +++-
src/getfem/getfem_modeling.h | 2 +-
src/getfem/getfem_models.h | 52 +-
src/getfem/getfem_nonlinear_elasticity.h | 9 +-
src/getfem/getfem_norm.h | 20 +-
src/getfem/getfem_plasticity.h | 12 +-
src/getfem_assembling_tensors.cc | 2 +-
src/getfem_contact_and_friction_common.cc | 8 +-
src/getfem_contact_and_friction_integral.cc | 59 +-
src/getfem_contact_and_friction_nodal.cc | 20 +-
src/getfem_generic_assembly.cc | 1155 ++++++++++++++++++++
src/getfem_level_set_contact.cc | 6 +-
src/getfem_mat_elem.cc | 32 +-
src/getfem_mesh.cc | 3 +-
src/getfem_mesh_fem_global_function.cc | 7 +-
src/getfem_mesh_region.cc | 6 +
src/getfem_model_solvers.cc | 3 +
src/getfem_models.cc | 1023 ++++++++++-------
src/gmm/gmm_MUMPS_interface.h | 29 +-
src/gmm/gmm_except.h | 54 +-
src/gmm/gmm_range_basis.h | 2 +-
src/gmm/gmm_sub_index.h | 29 +-
tests/elastostatic.cc | 8 +-
tests/elastostatic.param | 2 +-
tests/nonlinear_elastostatic.cc | 11 +-
68 files changed, 4811 insertions(+), 1025 deletions(-)
diff --cc configure.ac
index f02c135,7fd7a0f..920b17c
--- a/configure.ac
+++ b/configure.ac
@@@ -1265,8 -1194,6 +1295,14 @@@ case $host i
;;
esac
++<<<<<<< HEAD
+<<<<<<< HEAD:configure.in
+echo $shared_mode
+=======
+echo $shared_mode
+>>>>>>> upstream:configure.ac
++=======
+ if test "x$MSG" != "x"; then
+ echo -e "\n\nWARNINGS during the configure:\n$MSG\n\n"
+ fi
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
diff --cc interface/src/gf_cont_struct.cc
index 7f0e4f3,a31d4ea..1de9c3c
--- a/interface/src/gf_cont_struct.cc
+++ b/interface/src/gf_cont_struct.cc
@@@ -69,31 -52,6 +69,34 @@@ void gf_cont_struct(getfemint::mexargs_
(the default value is 'auto', which lets getfem choose itself);
possible values are 'superlu', 'mumps' (if supported), 'cg/ildlt',
'gmres/ilu' and 'gmres/ilut';
+<<<<<<< HEAD
++<<<<<<< HEAD
+ - 'max_iter', @int NIT
+ maximum number of iterations allowed in the correction (the default
+ value is 10);
+ - 'thr_iter', @int TIT
+ threshold number of iterations of the correction for enlarging the
+ step size (the default value is 8);
+ - 'max_res', @scalar RES
+ target residual value of the new point (the default value is 1e-6);
+ - 'max_diff', @scalar DIFF
+ determines a convergence criterion to the new tangent vector (the
+ default value is 1e-9);
+ - 'min_ang', @scalar ANG
+ minimal value of the cosine of the angle between tangents to the
+ solution curve at the old point and the new one (the default value
+ is 0.9);
+ - 'h_init', @scalar HIN
+ initial step size (the default value is 1e-2);
+ - 'h_max', @scalar HMAX
+ maximal step size (the default value is 1e-1);
+ - 'h_min', @scalar HMIN
+ minimal step size (the default value is 1e-5);
+=======
+ - 'bifurcations'
+ activates tools for detection and treatment of bifurcation points;
++=======
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
- 'h_init', @scalar HIN
initial step size (the default value is 1e-2);
- 'h_max', @scalar HMAX
diff --cc interface/src/scilab/demos/demo_continuation.sce
index 08c596c,3569c7d..4a6f241
--- a/interface/src/scilab/demos/demo_continuation.sce
+++ b/interface/src/scilab/demos/demo_continuation.sce
@@@ -1,9 -1,5 +1,13 @@@
// Scilab GetFEM++ interface
+<<<<<<< HEAD
++<<<<<<< HEAD
+// Copyright (C) 2011-2011 Tomas Ligursky, Yves Renard.
+=======
+// Copyright (C) 2011-2012 Tomas Ligursky, Yves Renard.
+>>>>>>> upstream
++=======
+ // Copyright (C) 2011-2013 Tomas Ligursky, Yves Renard.
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
//
// This file is a part of GetFEM++
//
@@@ -19,11 -15,8 +23,16 @@@
// along with this program; if not, write to the Free Software Foundation,
// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
+<<<<<<< HEAD
++<<<<<<< HEAD
+// Simple example of the bifurcation problem: -Delta(u) + u = lambda exp(u)
+=======
+// Simple example of the bifurcation problem: -Delta(u) + u = lambda * exp(u)
+>>>>>>> upstream
++=======
+ // Simple example of the bifurcation problem:
+ // -Delta(u) + u = lambda * exp(u).
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
//
// This program is used to check that scilab-getfem is working. This is also
// a good example of use of GetFEM++.
diff --cc interface/tests/matlab/demo_static_contact.m
index 95d22ab,e2cf479..0b05ec7
--- a/interface/tests/matlab/demo_static_contact.m
+++ b/interface/tests/matlab/demo_static_contact.m
@@@ -56,12 -52,11 +56,16 @@@ u_degree = 2
lambda_degree = 2;
incompressibility = 0;
p_degree = 1;
- r = 40; % Augmentation parameter
- gamma0 = 1/r; % Nitsche's method gamma0 parameter
- theta = 0; % Nitsche's method theta parameter
+ r = 40; % Augmentation parameter
+ gamma0 = 0.001; % Nitsche's method gamma0 parameter
+ theta = 1; % Nitsche's method theta parameter
++<<<<<<< HEAD
+condition_type = 0; % 0 = Explicitely kill horizontal rigid displacements
+>>>>>>> upstream
++=======
+ condition_type = 1; % 0 = Explicitely kill horizontal rigid displacements
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
% 1 = Kill rigid displacements using a global penalization
% 2 = Add a Dirichlet condition on the top of the structure
penalty_parameter = 1E-6; % Penalization coefficient for the global penalization
@@@ -75,11 -70,7 +79,15 @@@ end
niter = 100; % Maximum number of iterations for Newton's algorithm.
plot_mesh = true;
+<<<<<<< HEAD
++<<<<<<< HEAD
+version = 1; % 1 : frictionless contact and the basic contact brick
+=======
+version = 16; % 1 : frictionless contact and the basic contact brick
+>>>>>>> upstream
++=======
+ version = 16; % 1 : frictionless contact and the basic contact brick
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
% 2 : contact with 'static' Coulomb friction and basic contact brick
% 3 : frictionless contact and the contact with a
% rigid obstacle brick
diff --cc src/getfem/getfem_continuation.h
index 3f07ef2,069e4ad..1e97c7e
--- a/src/getfem/getfem_continuation.h
+++ b/src/getfem/getfem_continuation.h
@@@ -33,14 -33,7 +33,18 @@@
@author Tomas Ligursky <tomas.ligursky at gmail.com>
@author Yves Renard <Yves.Renard at insa-lyon.fr>
@date October 17, 2011.
++<<<<<<< HEAD
+ @brief (approximate) Moore-Penrose (also called Gauss-Newton) continuation method.
+<<<<<<< HEAD
+
+ NOTE: The bordered systems involved are solved by a block eliminiation
+ although the bordered matrix may be ill-conditioned in some cases!
+ Nevertheless, the algorithm seems to work well.
+=======
+>>>>>>> upstream
++=======
+ @brief inexact Moore-Penrose continuation method.
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
*/
#ifndef GETFEM_CONTINUATION_H__
#define GETFEM_CONTINUATION_H__
diff --cc src/getfem/getfem_models.h
index 8085b0d,c68af97..c4f7f80
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@@ -310,10 -279,7 +311,14 @@@ namespace getfem
const mimlist &mms, size_type reg)
: terms_to_be_computed(true), v_num(0), pbr(p), pdispatch(0), nbrhs(1),
vlist(vl), dlist(dl), tlist(tl), mims(mms), region(reg),
+<<<<<<< HEAD
++<<<<<<< HEAD
+ rveclist(1), rveclist_sym(1), cveclist(1),
+ cveclist_sym(1) { }
+=======
++=======
+ external_load(0),
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
partition( (mms.size()>0 ? &mms.at(0)->linked_mesh() : 0), region),
rveclist(1), rveclist_sym(1), cveclist(1),
cveclist_sym(1) { }
diff --cc src/getfem_contact_and_friction_integral.cc
index a70c143,0f41ead..f5ca05f
--- a/src/getfem_contact_and_friction_integral.cc
+++ b/src/getfem_contact_and_friction_integral.cc
@@@ -3467,240 -2643,47 +3467,271 @@@ namespace getfem
}
}
+ const bgeot::multi_index &sizes() const { return sizes_; }
- void compute(fem_interpolation_context &ctx, bgeot::base_tensor &t) {
+ virtual void compute(fem_interpolation_context&, bgeot::base_tensor &t);
+ virtual void prepare(fem_interpolation_context& /*ctx*/, size_type /*nb*/);
- md->compute_Neumann_terms(1, *varname, mf_u, U, ctx, n, tG);
- for (size_type i = 0; i < N; ++i)
- zeta[i] = tG[i]
- + ((g-un+alpha*un) * no[i] + alpha*wt[i] - alpha*u[i] ) / gamma;
- if ((option == 1) || (theta != scalar_type(0))) {
- coupled_projection(zeta, no, f_coeff, Pr);
- gmm::add(Pr, gmm::scaled(tG.as_vector(), -scalar_type(1)), pgg);
- }
+ };
- switch (option) {
- case 1:
- {
- ctx.pf()->real_base_value(ctx, tbv);
- size_type qmult = N / ctx.pf()->target_dim();
- size_type nbdofu = sizes_[0];
- if (theta != scalar_type(0)) {
- sizes_.resize(2);
- sizes_[1] = N;
- tp.adjust_sizes(sizes_);
- sizes_.resize(1);
- md->compute_Neumann_terms(2, *varname, mf_u, U, ctx, n, tp);
- }
- for (size_type i = 0; i < nbdofu; ++i) {
- t[i] = scalar_type(0);
- for (size_type j = 0; j < N; ++j) {
- if (theta != scalar_type(0))
- t[i] -= gamma*pgg[j]*theta*tp(i,j);
- if (qmult == 1) t[i] += Pr[j]*tbv(i,j);
- }
- if (qmult > 1) t[i] += Pr[i%N] * tbv(i/N,0);
- }
- }
- break;
++<<<<<<< HEAD
+ void contact_nitsche_nonlinear_term::compute
+ (fem_interpolation_context &/* ctx */, bgeot::base_tensor &t) {
- case 2:
- {
- size_type nbdofu = sizes_[1];
- sizes_[1] = N;
+ t.adjust_sizes(sizes_);
+ scalar_type e;
+ dim_type i, j, k, l;
+
+ if (option >= 3) { // computation of matrix A
+ e = f_coeff*gmm::neg(ln-r*(un-g));
+ auxN = lt - r*ut;
+ ball_projection_grad(auxN, e, GP);
+ ball_projection_grad_r(auxN, e, V);
+ e = Heav(r*(un-g) - ln);
+ gmm::rank_one_update(GP, no, gmm::scaled(V, -e*f_coeff));
+ gmm::rank_one_update(GP, gmm::scaled(no, e-gmm::vect_sp(GP,no,no)), no);
+ gmm::scale(GP, 1./r);
+ } else { // computation of vector W
+ e = gmm::neg(ln-r*(un-g));
+ V = lt - r*ut;
+ ball_projection(V, f_coeff*e);
+ V -= e*no;
+ }
+
+ switch (option) {
+ // one-dimensional tensors [N]
+ case 1:
+ for (i=0; i < N; ++i) t[i] = V[i];
+ break;
+
+ // two-dimensional tensors [N x N]
+ case 2:
+ V -= lnt;
+ gmm::scale(V, -1./r);
+ e = gmm::vect_sp(V, n);
+ for (i=0; i < N; ++i)
+ for (j=0; j < N; ++j) {
+ t(i,j) = mu*(V[i]*n[j]+V[j]*n[i]);
+ if (i == j) t(i,j) += lambda*e;
+ }
+ break;
+
+ case 3:
+ for (i=0; i < N; ++i)
+ for (j=0; j < N; ++j)
+ t(i,j) = r*r*GP(j,i);
+ break;
+
+ // three-dimensional tensors [N x N x N]
+ case 4:
+ gmm::mult(gmm::transposed(GP), n, V);
+ for (i=0; i < N; ++i)
+ for (j=0; j < N; ++j)
+ for (k=0; k < N; ++k) {
+ t(i,j,k) = -r*mu*(GP(j,i)*n[k] + GP(k,i)*n[j]);
+ if (j == k) t(i,j,k) -= r*lambda*V[i];
+ }
+ break;
+
+ case 5:
+ gmm::mult(GP, n, V);
+ for (i=0; i < N; ++i)
+ for (j=0; j < N; ++j)
+ for (k=0; k < N; ++k) {
+ t(i,j,k) = -r*mu*(GP(k,i)*n[j] + GP(k,j)*n[i]);
+ if (i == j) t(i,j,k) -= r*lambda*V[k];
+ }
+ break;
+
+ // four-dimensional tensors [N x N x N x N]
+
+ case 6:
+
+ for (i=0; i < N; ++i) GP(i,i) -= 1./r; // matrix B
+
+ e = gmm::vect_sp(GP, n, n);
+ gmm::mult(gmm::transposed(GP), n, auxN);
+ gmm::mult(GP, n, V);
+
+ for (i=0; i < N; ++i)
+ for (j=0; j < N; ++j)
+ for (k=0; k < N; ++k)
+ for (l=0; l < N; ++l) {
+ t(i,j,k,l) = mu*mu*(n[i]*GP(k,j)*n[l] + n[j]*GP(k,i)*n[l]
+ + n[j]*GP(l,i)*n[k] + n[i]*GP(l,j)*n[k]);
+ if (i == j && k == l) t(i,j,k,l) += lambda*lambda*e;
+ if (i == j) t(i,j,k,l) += lambda*mu*(V[k]*n[l] + V[l]*n[k]);
+ if (k == l) t(i,j,k,l) += lambda*mu*(auxN[j]*n[i]+auxN[i]*n[j]);
+ }
+
+ break;
+ default : GMM_ASSERT1(false, "Invalid option");
+ }
+ }
+
+
+ void contact_nitsche_nonlinear_term::prepare
+ (fem_interpolation_context& ctx, size_type nb) {
+ size_type cv = ctx.convex_num();
+
+ switch (nb) { // last is computed first
+ case 1 : // calculate [un] and [ut] interpolating [U] on [mf_u]
+ coeff.resize(mf_u.nb_basic_dof_of_element(cv));
+ gmm::copy(gmm::sub_vector(U, gmm::sub_index
+ (mf_u.ind_basic_dof_of_element(cv))), coeff);
+ ctx.pf()->interpolation(ctx, coeff, V, N);
+ un = gmm::vect_sp(V, no);
+ ut = V - un * no;
+ ctx.pf()->interpolation_grad(ctx, coeff, GP, N);
+ lnt = lambda*(gmm::mat_trace(GP))*n;
+ gmm::mult_add(GP, gmm::scaled(n, mu), lnt);
+ gmm::mult_add(gmm::transposed(GP), gmm::scaled(n, mu), lnt);
+ ln = gmm::vect_sp(lnt, no);
+ lt = lnt - ln * no;
+ break;
+
+ case 2 : // calculate [g] and [no] interpolating [obs] on [mf_obs]
+ // calculate [ln] and [lt] from [lnt] and [no]
+ coeff.resize(mf_obs.nb_basic_dof_of_element(cv));
+ gmm::copy(gmm::sub_vector
+ (obs, gmm::sub_index
+ (mf_obs.ind_basic_dof_of_element(cv))), coeff);
+ ctx.pf()->interpolation_grad(ctx, coeff, grad, 1);
+ gmm::copy(gmm::mat_row(grad, 0), no);
+ no /= -gmm::vect_norm2(no);
+ ctx.pf()->interpolation(ctx, coeff, aux1, 1);
+ g = aux1[0];
+ n = bgeot::compute_normal(ctx, ctx.face_num());
+ n /= gmm::vect_norm2(n);
+ break;
+
+ case 3 :// calculate [f_coeff] interpolating [friction_coeff] on [mf_coeff]
+ if (pmf_coeff) {
+ coeff.resize(pmf_coeff->nb_basic_dof_of_element(cv));
+ gmm::copy(gmm::sub_vector
+ (friction_coeff, gmm::sub_index
+ (pmf_coeff->ind_basic_dof_of_element(cv))), coeff);
+ ctx.pf()->interpolation(ctx, coeff, aux1, 1);
+ f_coeff = aux1[0];
+ }
+ break;
+
+ default : GMM_ASSERT1(false, "Invalid option");
+ }
+ }
+
+
+
+
+ template<typename MAT, typename VECT1>
+ void asm_Nitsche_contact_rigid_obstacle_tangent_matrix
+ (MAT &K, const mesh_im &mim,
+ const getfem::mesh_fem &mf_u, const VECT1 &U,
+ const getfem::mesh_fem &mf_obs, const VECT1 &obs,
+ const getfem::mesh_fem *pmf_coeff, const VECT1 &f_coeff,
+ scalar_type gamma, scalar_type lambda, scalar_type mu,
+ const mesh_region &rg, int option = 1) {
+
+ contact_nitsche_nonlinear_term
+ nterm1(6, gamma, lambda, mu, mf_u, U, mf_obs, obs, pmf_coeff, &f_coeff),
+ nterm2(3, gamma, lambda, mu, mf_u, U, mf_obs, obs, pmf_coeff, &f_coeff),
+ nterm3(4, gamma, lambda, mu, mf_u, U, mf_obs, obs, pmf_coeff, &f_coeff),
+ nterm4(5, gamma, lambda, mu, mf_u, U, mf_obs, obs, pmf_coeff, &f_coeff);
+
+ const std::string aux_fems = pmf_coeff ? "#1,#2,#3" : "#1,#2";
+
+ getfem::generic_assembly assem;
+ std::string as_str
+ = ((option == 0) ? "w1=comp(NonLin$1(#1,"+aux_fems+")(i,j,k,l).vGrad(#1)(:,i,j).vGrad(#1)(:,k,l));" : "")
+ + "w2=comp(NonLin$2(#1,"+aux_fems+").vBase(#1).vBase(#1))(i,j,:,i,:,j);"
+ + "w3=comp(NonLin$3(#1,"+aux_fems+").vBase(#1).vGrad(#1))(i,j,k,:,i,:,j,k);"
+ + ((option == 0) ? "w4=comp(NonLin$4(#1,"+aux_fems+").vGrad(#1).vBase(#1))(i,j,k,:,i,j,:,k);" : "")
+ + ((option == 0) ? "M(#1,#1)+=w1+w2+w3+w4;" : "M(#1,#1)+=w2+w3;");
+
+ assem.set(as_str);
+=======
+ if (f_coeff_) f_coeff = (*f_coeff_)[0]; else f_coeff = scalar_type(0);
+ else {
+ friction_coeff.resize(pmf_coeff->nb_basic_dof());
+ pmf_coeff->extend_vector(*f_coeff_, friction_coeff);
+ }
+ if (WT_) {
+ WT.resize(mf_u.nb_basic_dof());
+ mf_u_.extend_vector(*WT_, WT);
+ }
+ }
+
+
+ void compute(fem_interpolation_context &ctx, bgeot::base_tensor &t) {
+
+ md->compute_Neumann_terms(1, *varname, mf_u, U, ctx, n, tG);
+ for (size_type i = 0; i < N; ++i)
+ zeta[i] = tG[i]
+ + ((g-un+alpha*un) * no[i] + alpha*wt[i] - alpha*u[i] ) / gamma;
+ if ((option == 1) || (theta != scalar_type(0))) {
+ coupled_projection(zeta, no, f_coeff, Pr);
+ gmm::add(Pr, gmm::scaled(tG.as_vector(), -scalar_type(1)), pgg);
+ }
+
+ switch (option) {
+ case 1:
+ {
+ ctx.pf()->real_base_value(ctx, tbv);
+ size_type qmult = N / ctx.pf()->target_dim();
+ short_type nbdofu = sizes_[0];
+ if (theta != scalar_type(0)) {
+ sizes_.resize(2);
+ sizes_[1] = N;
+ tp.adjust_sizes(sizes_);
+ sizes_.resize(1);
+ md->compute_Neumann_terms(2, *varname, mf_u, U, ctx, n, tp);
+ }
+ for (size_type i = 0; i < nbdofu; ++i) {
+ t[i] = scalar_type(0);
+ for (size_type j = 0; j < N; ++j) {
+ if (theta != scalar_type(0))
+ t[i] -= gamma*pgg[j]*theta*tp(i,j);
+ if (qmult == 1) t[i] += Pr[j]*tbv(i,j);
+ }
+ if (qmult > 1) t[i] += Pr[i%N] * tbv(i/N,0);
+ }
+ }
+ break;
+
+ case 2:
+ {
+ short_type nbdofu = sizes_[1];
++=======
++ switch (option) {
++ case 1:
++ {
++ ctx.pf()->real_base_value(ctx, tbv);
++ size_type qmult = N / ctx.pf()->target_dim();
++ size_type nbdofu = sizes_[0];
++ if (theta != scalar_type(0)) {
++ sizes_.resize(2);
++ sizes_[1] = N;
++ tp.adjust_sizes(sizes_);
++ sizes_.resize(1);
++ md->compute_Neumann_terms(2, *varname, mf_u, U, ctx, n, tp);
++ }
++ for (size_type i = 0; i < nbdofu; ++i) {
++ t[i] = scalar_type(0);
++ for (size_type j = 0; j < N; ++j) {
++ if (theta != scalar_type(0))
++ t[i] -= gamma*pgg[j]*theta*tp(i,j);
++ if (qmult == 1) t[i] += Pr[j]*tbv(i,j);
++ }
++ if (qmult > 1) t[i] += Pr[i%N] * tbv(i/N,0);
++ }
++ }
++ break;
++
++ case 2:
++ {
++ size_type nbdofu = sizes_[1];
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
+ sizes_[1] = N;
tp.adjust_sizes(sizes_);
sizes_[1] = nbdofu;
md->compute_Neumann_terms(2, *varname, mf_u, U, ctx, n, tp);
diff --cc src/getfem_contact_and_friction_nodal.cc
index e64c9bb,21d3caa..6d84157
--- a/src/getfem_contact_and_friction_nodal.cc
+++ b/src/getfem_contact_and_friction_nodal.cc
@@@ -1208,11 -1173,6 +1212,14 @@@ namespace getfem
bool Hughes_stabilized_=false,
bool friction_dynamic_term_=false) {
++<<<<<<< HEAD
+#if GETFEM_PARA_LEVEL > 1
+ if (!getfem::MPI_IS_MASTER()) GMM_WARNING1("Nodal contact bricks don't support GETFEM_PARA_LEVEL > 1 yet!!!");
+#endif
+
+>>>>>>> upstream:src/getfem_contact_and_friction_nodal.cc
++=======
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
if (aug_version == 4 && contact_only_) aug_version = 3;
augmentation_version = aug_version;
GMM_ASSERT1(aug_version >= 1 && aug_version <= 4,
diff --cc src/getfem_model_solvers.cc
index 313a148,7ce7ac6..d71f073
--- a/src/getfem_model_solvers.cc
+++ b/src/getfem_model_solvers.cc
@@@ -85,312 -85,6 +85,315 @@@ namespace getfem
}
/* ***************************************************************** */
++<<<<<<< HEAD
+ /* Intermediary structure for Newton algorithms. */
+ /* ***************************************************************** */
+
+ #define TRACE_SOL 0
+
+ template <typename MAT, typename VEC>
+ struct model_pb {
+
+ typedef MAT MATRIX;
+ typedef VEC VECTOR;
+ typedef typename gmm::linalg_traits<VECTOR>::value_type T;
+ typedef typename gmm::number_traits<T>::magnitude_type R;
+
+ model &md;
+ bool is_reduced;
+ std::vector<size_type> &sind;
+ gmm::sub_index I;
+ abstract_newton_line_search &ls;
+ VECTOR stateinit, &state;
+ const VECTOR &rhs;
+ const MATRIX &K;
+ MATRIX Kr;
+ VECTOR rhsr;
+ bool with_pseudo_potential;
+
+ void compute_tangent_matrix(void) {
+ md.to_variables(state);
+ md.assembly(model::BUILD_MATRIX);
+ if (is_reduced) {
+ gmm::resize(Kr, sind.size(), sind.size());
+ gmm::copy(gmm::sub_matrix(K, I, I), Kr);
+ }
+ }
+
+ const MATRIX &tangent_matrix(void) { return (is_reduced ? Kr : K); }
+
+ inline T scale_residual(void) const { return T(1); }
+
+ void compute_residual(void) {
+ md.to_variables(state);
+ md.assembly(model::BUILD_RHS);
+ if (is_reduced) {
+ gmm::resize(rhsr, sind.size());
+ gmm::copy(gmm::sub_vector(rhs, I), rhsr);
+ }
+ }
+
+ void compute_pseudo_potential(void)
+ { md.to_variables(state); md.assembly(model::BUILD_PSEUDO_POTENTIAL); }
+
+ void perturbation(void) {
+ R res = gmm::vect_norm2(state), ampl = std::max(res*R(1E-20), R(1E-50));
+ std::vector<R> V(gmm::vect_size(state));
+ gmm::fill_random(V);
+ gmm::add(gmm::scaled(V, ampl), state);
+ }
+
+ const VECTOR &residual(void) { return (is_reduced ? rhsr : rhs); }
+
+ R residual_norm(void) { // A norm1 seems to be better than a norm2
+ // at least for contact problems.
+ if (is_reduced)
+ return gmm::vect_norm1(rhsr)/R(gmm::vect_size(rhsr));
+ else
+ return gmm::vect_norm1(rhs)/R(gmm::vect_size(rhs));
+ }
+
+ R compute_res(bool comp = true) {
+ if (with_pseudo_potential) {
+ compute_pseudo_potential();
+ return md.pseudo_potential();
+ } else {
+ if (comp) compute_residual();
+ return residual_norm();
+<<<<<<< HEAD
+ }
+ }
+
+
+ R line_search(VECTOR &dr, const gmm::iteration &iter) {
+ size_type nit = 0;
+ gmm::resize(stateinit, md.nb_dof());
+ gmm::copy(state, stateinit);
+ R alpha(1), res, res_init, R0;
+
+ res_init = res = compute_res(false);
+ // cout << "first residual = " << residual() << endl << endl;
+ R0 = (is_reduced ? gmm::real(gmm::vect_sp(dr, rhsr))
+ : gmm::real(gmm::vect_sp(dr, rhs)));
+
+ // store the initial residual for the reprojection step.
+// size_type N = gmm::vect_size(residual());
+// VECTOR first_rhs(N); gmm::copy(residual(), first_rhs); // � �liminer
+
+
+ // Compute the second derivative at alpha = 0 (value 2*a)
+ // not very effective ... precision problem ?
+
+// R EPS = 1e-6;
+// gmm::add(gmm::sub_vector(stateinit, I), gmm::scaled(dr, EPS),
+// gmm::sub_vector(state, I));
+// R res2 = compute_res();
+
+// R a = (res2 + EPS * res - res) / gmm::sqr(EPS);
+// R delta = gmm::sqr(res)-R(4)*a*res;
+// R alpha_bis = R(1);
+// if (delta > 0) {
+// R l1 = (res - sqrt(delta)) / (R(2)*a);
+// R l2 = (res - sqrt(delta)) / (R(2)*a);
+// if (l1 > 0) alpha_bis = std::min(alpha_bis, l1);
+// if (l2 > 0) alpha_bis = std::min(alpha_bis, l2);
+// }
+// if (a > 0)
+// alpha_bis = std::min(alpha_bis, res / (R(2)*a));
+// cout << "alpha_bis = " << alpha_bis << endl;
+
+// gmm::add(gmm::sub_vector(stateinit, I), gmm::scaled(dr, alpha_bis),
+// gmm::sub_vector(state, I));
+
+// R res3 = compute_res();
+// cout << "res3 = " << res3 << endl;
+
+// cout.precision(16);
+// cout << "res = " << res << " res2 = " << res2 << endl;
+// cout << "a = " << (res2 + EPS * res - res) / gmm::sqr(EPS) << endl;
+// cout << "aa = " << (res2 + EPS * res - res) << endl;
+
+
+=======
+ }
+ }
+
+
+ R line_search(VECTOR &dr, const gmm::iteration &iter) {
+ size_type nit = 0;
+ gmm::resize(stateinit, md.nb_dof());
+ gmm::copy(state, stateinit);
+ R alpha(1), res, /* res_init, */ R0;
+
+ /* res_init = */ res = compute_res(false);
+ // cout << "first residual = " << residual() << endl << endl;
+ R0 = (is_reduced ? gmm::real(gmm::vect_sp(dr, rhsr))
+ : gmm::real(gmm::vect_sp(dr, rhs)));
+
+ // store the initial residual for the reprojection step.
+// size_type N = gmm::vect_size(residual());
+// VECTOR first_rhs(N); gmm::copy(residual(), first_rhs); // � �liminer
+
+
+ // Compute the second derivative at alpha = 0 (value 2*a)
+ // not very effective ... precision problem ?
+
+// R EPS = 1e-6;
+// gmm::add(gmm::sub_vector(stateinit, I), gmm::scaled(dr, EPS),
+// gmm::sub_vector(state, I));
+// R res2 = compute_res();
+
+// R a = (res2 + EPS * res - res) / gmm::sqr(EPS);
+// R delta = gmm::sqr(res)-R(4)*a*res;
+// R alpha_bis = R(1);
+// if (delta > 0) {
+// R l1 = (res - sqrt(delta)) / (R(2)*a);
+// R l2 = (res - sqrt(delta)) / (R(2)*a);
+// if (l1 > 0) alpha_bis = std::min(alpha_bis, l1);
+// if (l2 > 0) alpha_bis = std::min(alpha_bis, l2);
+// }
+// if (a > 0)
+// alpha_bis = std::min(alpha_bis, res / (R(2)*a));
+// cout << "alpha_bis = " << alpha_bis << endl;
+
+// gmm::add(gmm::sub_vector(stateinit, I), gmm::scaled(dr, alpha_bis),
+// gmm::sub_vector(state, I));
+
+// R res3 = compute_res();
+// cout << "res3 = " << res3 << endl;
+
+// cout.precision(16);
+// cout << "res = " << res << " res2 = " << res2 << endl;
+// cout << "a = " << (res2 + EPS * res - res) / gmm::sqr(EPS) << endl;
+// cout << "aa = " << (res2 + EPS * res - res) << endl;
+
+#if TRACE_SOL
+ static int trace_number = 0;
+ int trace_iter = 0;
+ {
+ std::stringstream trace_name;
+ trace_name << "line_search_state" << std::setfill('0')
+ << std::setw(3) << trace_number << "_000_init";
+ gmm::vecsave(trace_name.str(),stateinit);
+ }
+ trace_number++;
+#endif
+
+>>>>>>> upstream
+ ls.init_search(res, iter.get_iteration(), R0);
+ do {
+ alpha = ls.next_try();
+ gmm::add(gmm::sub_vector(stateinit, I), gmm::scaled(dr, alpha),
+ gmm::sub_vector(state, I));
+<<<<<<< HEAD
+=======
+#if TRACE_SOL
+ {
+ trace_iter++;
+ std::stringstream trace_name;
+ trace_name << "line_search_state" << std::setfill('0')
+ << std::setw(3) << trace_number << "_"
+ << std::setfill('0') << std::setw(3) << trace_iter;
+ gmm::vecsave(trace_name.str(), state);
+ }
+#endif
+>>>>>>> upstream
+ res = compute_res();
+ // cout << "residual = " << residual() << endl << endl;
+ R0 = (is_reduced ? gmm::real(gmm::vect_sp(dr, rhsr))
+ : gmm::real(gmm::vect_sp(dr, rhs)));
+
+ ++ nit;
+ } while (!ls.is_converged(res, R0));
+
+ if (alpha != ls.converged_value() || with_pseudo_potential) {
+ alpha = ls.converged_value();
+ gmm::add(gmm::sub_vector(stateinit, I), gmm::scaled(dr, alpha),
+ gmm::sub_vector(state, I));
+ res = ls.converged_residual();
+ compute_residual();
+ }
+
+// if (1) { // reprojection step ... only for non reduced case ...
+// GMM_ASSERT1(!is_reduced, "to be done ...");
+// // detection des ddls "coupables"
+// R mean = R(0);
+// for (size_type i = 0; i < N; ++i) mean += gmm::abs(residual()[i]);
+// mean /= R(N);
+// cout << "mean = " << mean << endl << "selected : ";
+
+// VECTOR new_dr(N);
+// size_type Ndof = 0;
+// for (size_type i = 0; i < N; ++i)
+// if (gmm::abs(residual()[i])
+// > std::max(gmm::abs(first_rhs[i]), mean) * R(4)) {
+// cout << i << " ";
+// new_dr[i] = dr[i];
+// ++Ndof;
+// }
+// cout << endl;
+
+// if (Ndof > 0) {
+
+// R EPS = 1e-8;
+// R THR = 0.8;
+// R mini_alpha = alpha * THR;
+
+// cout << "performing post-correction" << endl;
+// cout << "res = " << res << " res_init = " << res_init << endl;
+
+
+// // divided difference
+// gmm::add(gmm::scaled(new_dr, -EPS), state);
+// R res2 = compute_res();
+// cout << "pente = " << (res - res2) / EPS << endl;
+// R beta = (res - res2) / EPS;
+// mini_alpha = (res2 - THR*std::min(res_init, res))/((beta != R(0)) ? beta : R(1));
+
+// cout << "mini_alpha = " << mini_alpha << endl;
+// if (mini_alpha > alpha * THR)
+// mini_alpha = alpha * THR;
+
+// if (mini_alpha > 0) {
+
+// gmm::add(gmm::scaled(new_dr, -mini_alpha+EPS), state);
+// R res4 = compute_res(), res3;
+// cout << "res4 = " << res4 << endl;
+// size_type mini_it = 0;
+// do {
+// res3 = res4;
+// mini_alpha /= R(2); mini_it++;
+// gmm::add(gmm::scaled(new_dr, mini_alpha), state);
+// res4 = compute_res();
+// cout << "res4 = " << res4 << endl;
+// } while (res4 < res3 && mini_it <= 6);
+
+// if (res4 < res3) { res3 = res4; }
+// else gmm::add(gmm::scaled(new_dr, -mini_alpha), state);
+
+// if (res3 > res) { // cancel correction
+// gmm::add(stateinit, gmm::scaled(dr, alpha), state);
+// }
+// }
+// else gmm::add(stateinit, gmm::scaled(dr, alpha), state);
+// res = compute_res();
+// }
+// }
+ return alpha;
+ }
+
+ model_pb(model &m, abstract_newton_line_search &ls_, VECTOR &st,
+ const VECTOR &rhs_, const MATRIX &K_, bool reduced_,
+ std::vector<size_type> &sind_,
+ bool with_pseudo_pot = false)
+ : md(m), is_reduced(reduced_), sind(sind_), I(sind_), ls(ls_), state(st),
+ rhs(rhs_), K(K_), with_pseudo_potential(with_pseudo_pot) {}
+
+ };
+
+ /* ***************************************************************** */
++=======
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
/* Standard solve. */
/* ***************************************************************** */
diff --cc src/getfem_models.cc
index e635d9a,d607e94..7cec4a4
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@@ -1,10 -1,7 +1,10 @@@
+<<<<<<< HEAD
+=======
/* -*- c++ -*- (enables emacs c++ mode) */
+>>>>>>> upstream
/*===========================================================================
- Copyright (C) 2009-2012 Yves Renard
+ Copyright (C) 2009-2013 Yves Renard
This file is a part of GETFEM++
@@@ -213,10 -215,8 +224,13 @@@ namespace getfem
// the primal variable. A search is done on all the terms of the
// model. Only the the corresponding linear terms are added.
// If no linear term is available, a mass matrix is used.
++<<<<<<< HEAD
+>>>>>>> upstream
+
++=======
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
gmm::col_matrix< gmm::rsvector<scalar_type> >
- MM(it2->second.mf->nb_dof(), it->second.mf->nb_dof());
+ MM(it2->second.associated_mf().nb_dof(), it->second.mf->nb_dof());
bool termadded = false;
if (it->second.filter == VDESCRFILTER_CTERM) {
@@@ -272,26 -263,15 +286,34 @@@
if (!termadded)
GMM_WARNING1("No term found to filter multiplier " << it->first
<< ". Variable is cancelled");
++<<<<<<< HEAD
+#if GETFEM_PARA_LEVEL > 1
+ if (termadded) {
+ // we assume that all bricks take mpi_region into account but it
+ // would be better if the brick itself could report if it supports
+ // distributed assembly
+ // This is only a reference implementation, it needs to be optimized
+ // maybe by using gmm::mpi_distributed_matrix
+ std::vector<scalar_type> tmpvec1(gmm::mat_nrows(MM)), tmpvec2(gmm::mat_nrows(MM));
+ for (size_type k = 0; k < gmm::mat_ncols(MM); ++k) {
+ gmm::copy(gmm::mat_col(MM,k),tmpvec1);
+ MPI_SUM_VECTOR(tmpvec1,tmpvec2);
+ gmm::copy(tmpvec2,gmm::mat_col(MM,k));
+ }
+ }
+#endif
+>>>>>>> upstream
++=======
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
} else if (it->second.filter == VDESCRFILTER_INFSUP) {
- asm_mass_matrix(MM, *(it->second.mim), *(it2->second.mf),
- *(it->second.mf), it->second.m_region);
+ mesh_region rg(it->second.m_region);
+ it->second.mim->linked_mesh().intersect_with_mpi_region(rg);
+ asm_mass_matrix(MM, *(it->second.mim), it2->second.associated_mf(),
+ *(it->second.mf), rg);
}
+
+ MPI_SUM_SPARSE_MATRIX(MM);
+
//
// filtering
//
@@@ -316,13 -297,20 +339,24 @@@
}
}
+ // #if GETFEM_PARA_LEVEL > 1
+ // if (!rk) cout << "Range basis for multipliers for " << itbd->first << " time : " << MPI_Wtime()-tt_ref << endl;
+
+ // #endif
+
if (mults.size() > 1) {
+<<<<<<< HEAD
+ range_basis(MGLOB, glob_columns, 1E-12, gmm::col_major(), true);
+=======
gmm::range_basis(MGLOB, glob_columns, 1E-12, gmm::col_major(), true);
+>>>>>>> upstream
+
+ // #if GETFEM_PARA_LEVEL > 1
+ // if (!rk) cout << "Producing partial mf for multipliers for " << itbd->first << " time : " << MPI_Wtime()-tt_ref << endl;
+
+ // #endif
+
s = 0;
for (size_type k = 0; k < mults.size(); ++k) {
VAR_SET::iterator it = variables.find(mults[k]);
@@@ -1661,9 -1502,9 +1711,11 @@@
// brick.rmatlist = real_matlist(brick.tlist.size());
// brick.rveclist[0] = real_veclist(brick.tlist.size());
// }
+<<<<<<< HEAD
+=======
+
+ if (version & BUILD_RHS) approx_external_load_ += brick.external_load;
}
if (version & BUILD_RHS) {
@@@ -1803,16 -1666,14 +1878,18 @@@
const model_real_plain_vector &
model::real_variable(const std::string &name, size_type niter) const {
GMM_ASSERT1(!complex_version, "This model is a complex one");
- context_check(); if (act_size_to_be_done) actualize_sizes();
+ context_check();
VAR_SET::const_iterator it = variables.find(name);
GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name);
+ if (act_size_to_be_done && it->second.filter != VDESCRFILTER_NO)
+ actualize_sizes();
if (niter == size_type(-1)) niter = it->second.default_iter;
GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
+<<<<<<< HEAD
+ "Unvalid iteration number "
+=======
"Invalid iteration number "
+>>>>>>> upstream
<< niter << " for " << name);
return it->second.real_value[niter];
}
@@@ -1820,16 -1681,14 +1897,18 @@@
const model_complex_plain_vector &
model::complex_variable(const std::string &name, size_type niter) const {
GMM_ASSERT1(complex_version, "This model is a real one");
- context_check(); if (act_size_to_be_done) actualize_sizes();
+ context_check();
VAR_SET::const_iterator it = variables.find(name);
GMM_ASSERT1(it!=variables.end(), "Undefined variable " << name);
+ if (act_size_to_be_done && it->second.filter != VDESCRFILTER_NO)
+ actualize_sizes();
if (niter == size_type(-1)) niter = it->second.default_iter;
GMM_ASSERT1(it->second.n_iter + it->second.n_temp_iter > niter,
+<<<<<<< HEAD
+ "Unvalid iteration number "
+=======
"Invalid iteration number "
+>>>>>>> upstream
<< niter << " for " << name);
return it->second.complex_value[niter];
}
@@@ -1869,20 -1724,15 +1952,32 @@@
<< niter << " for " << name);
return it->second.complex_value[niter];
}
++<<<<<<< HEAD
+ void model::check_brick_stiffness_rhs(size_type ind_brick) const
+ {
+
+
+ const brick_description &brick = bricks[ind_brick];
+ update_brick(ind_brick, model::BUILD_ALL);
+
+<<<<<<< HEAD
+ brick.pbr->check_stiffness_matrix_and_rhs(*this, ind_brick,
+=======
+ brick.pbr->check_stiffness_matrix_and_rhs(*this, ind_brick, brick.tlist,
+>>>>>>> upstream
+ brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
+ brick.rveclist[0], brick.rveclist_sym[0], brick.region);
++=======
+
+ void model::check_brick_stiffness_rhs(size_type ind_brick) const {
+
+ const brick_description &brick = bricks[ind_brick];
+ update_brick(ind_brick, model::BUILD_ALL);
+
+ brick.pbr->check_stiffness_matrix_and_rhs(*this, ind_brick, brick.tlist,
+ brick.vlist, brick.dlist, brick.mims, brick.rmatlist,
+ brick.rveclist[0], brick.rveclist_sym[0], brick.region);
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
}
@@@ -1893,222 -1743,139 +1988,318 @@@
//
//
// ----------------------------------------------------------------------
+<<<<<<< HEAD
++<<<<<<< HEAD
+ void virtual_brick::check_stiffness_matrix_and_rhs(const model &md, size_type s,
+ const model::varnamelist &vl,
+ const model::varnamelist &dl,
+ const model::mimlist &mims,
+ model::real_matlist &matl,
+ model::real_veclist &rvc1,
+ model::real_veclist &rvc2,
+ size_type rg,
+ const scalar_type TINY) const
+ {
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2, rg, model::BUILD_MATRIX);
+ model_real_sparse_matrix SM(matl[0]);
+ gmm::fill(rvc1[0], 0.0);
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2, rg, model::BUILD_RHS);
+ model_real_plain_vector RHS0(rvc1[0]);
+
+ //finite difference stiffness
+ model_real_sparse_matrix fdSM(matl[0].nrows(),matl[0].ncols());
+
+ for (size_type i=0;i<rvc1[0].size();i++){
+ model_real_plain_vector U(md.real_variable(vl[0]));
+ U[i]+=TINY;
+ gmm::copy(U, md.set_real_variable(vl[0]));
+ gmm::fill(rvc1[0], 0.0);
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2, rg, model::BUILD_RHS);
+ model_real_plain_vector RHS1(rvc1[0]);
+ for (size_type j=0;j<rvc1[0].size();j++){
+ fdSM(i,j)=(RHS0[j]-RHS1[j])/TINY;
+ }
+ U[i]-=TINY;
+ gmm::copy(U, md.set_real_variable(vl[0]));
+ }
+ model_real_sparse_matrix diffSM(matl[0].nrows(),matl[0].ncols());
+ gmm::add(matl[0],gmm::scaled(fdSM,-1.0),diffSM);
+ scalar_type norm_error_euc = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(matl[0])*100;
+ scalar_type norm_error_1 = gmm::mat_norm1(diffSM)/gmm::mat_norm1(matl[0])*100;
+ scalar_type norm_error_max = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(matl[0])*100;
+
+ model_real_sparse_matrix diffSMtransposed(matl[0].nrows(),matl[0].ncols());
+ gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
+ scalar_type nsym_norm_error_euc = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
+ scalar_type nsym_norm_error_1 = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
+ scalar_type nsym_norm_error_max = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
+
+ //print matrix if the size is small
+ if(rvc1[0].size()<8){
+ std::cout << "RHS Stiffness Matrix: \n";
+ std::cout << "------------------------\n";
+ for(size_type i=0; i < rvc1[0].size(); ++i){
+ std::cout << "[";
+ for(size_type j=0; j < rvc1[0].size(); ++j){
+ std::cout << fdSM(i,j) << " ";
+ }
+ std::cout << "]\n";
+ }
+ std::cout << "Analytical Stiffness Matrix: \n";
+ std::cout << "------------------------\n";
+ for(size_type i=0; i < rvc1[0].size(); ++i){
+ std::cout << "[";
+ for(size_type j=0; j < rvc1[0].size(); ++j){
+ std::cout << matl[0](i,j) << " ";
+ }
+ std::cout << "]\n";
+ }
+ std::cout << "Vector U: \n";
+ std::cout << "------------------------\n";
+ for(size_type i=0; i < rvc1[0].size(); ++i){
+ std::cout << "[";
+ std::cout << md.real_variable(vl[0])[i] << " ";
+ std::cout << "]\n";
+ }
+ }
+
+ std::cout<<"\n\nfinite diff test error_norm_eucl: "<<norm_error_euc <<"%"<<std::endl;
+ std::cout<<"finite diff test error_norm1: "<<norm_error_1 <<"%"<<std::endl;
+ std::cout<<"finite diff test error_max_norm: "<<norm_error_max <<"%"<<std::endl;
+ std::cout<<"\n\nNonsymmetrical test error_norm_eucl: "<<nsym_norm_error_euc <<"%"<<std::endl;
+ std::cout<<"Nonsymmetrical test error_norm1: "<<nsym_norm_error_1 <<"%"<<std::endl;
+ std::cout<<"Nonsymmetrical test error_max_norm: "<<nsym_norm_error_max <<"%"<<std::endl;
+ }
+
+=======
+ void virtual_brick::check_stiffness_matrix_and_rhs
+ (const model &md, size_type s,
+ const model::termlist& tlist,
+ const model::varnamelist &vl,
+ const model::varnamelist &dl,
+ const model::mimlist &mims,
+ model::real_matlist &matl,
+ model::real_veclist &rvc1,
+ model::real_veclist &rvc2,
+ size_type rg,
+ const scalar_type TINY) const {
+ std::cout<<"******Verifying stiffnesses of *******"<<std::endl;
+ std::cout<<"*** "<<brick_name()<<std::endl;
+
+ //Build the index for the corresponding RHS
+ std::map<std::string,size_type> rhs_index;
+ for(size_type iterm=0;iterm<matl.size();iterm++)
+ if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
+
+ if (rhs_index.size()==0){
+ GMM_WARNING0("*** cannot verify stiffness for this brick***");
+ return;
+ }
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2,
+ rg, model::BUILD_MATRIX);
+ for(size_type iterm=0;iterm<matl.size();iterm++){
+
+ std::cout<<std::endl;
+ std::cout<<" Stiffness["<<tlist[iterm].var1
+ <<","<<tlist[iterm].var2<<"]:"<<std::endl;
+ if (md.real_variable(tlist[iterm].var1).size()==0)
+ {
+ std::cout<<" "<<tlist[iterm].var1<<" has zero size. Skipping this term"<<std::endl;
+ continue;
+ }
+ if (md.real_variable(tlist[iterm].var2).size()==0)
+ {
+ std::cout<<" "<<tlist[iterm].var2<<" has zero size. Skipping this term"<<std::endl;
+ continue;
+ }
+
+ model_real_sparse_matrix SM(matl[iterm]);
+ gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2,
+ rg, model::BUILD_RHS);
+ if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
+ std::cout<<" The assembled matrix is nearly zero, skipping."<<std::endl;
+ continue;
+ }
+ model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
+
+ //finite difference stiffness
+ model_real_sparse_matrix fdSM(matl[iterm].nrows(),matl[iterm].ncols());
+ model_real_plain_vector&U = md.set_real_variable(tlist[iterm].var2);
+ model_real_plain_vector& RHS1 =rvc1[rhs_index[tlist[iterm].var1]];
+ for (size_type j=0; j < matl[iterm].ncols(); j++){
+ U[j]+=TINY;
+ gmm::fill(RHS1, 0.0);
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2,
+ rg, model::BUILD_RHS);
+ for (size_type i=0;i<matl[iterm].nrows();i++)
+ fdSM(i,j) = (RHS0[i]-RHS1[i])/TINY;
+ U[j]-=TINY;
+ }
+
+ model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
+ gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
+ scalar_type norm_error_euc
+ = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
+ scalar_type norm_error_1
+ = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
+ scalar_type norm_error_max
+ = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
+
+ //checking symmetry of diagonal terms
+ scalar_type nsym_norm_error_euc=0.0;
+ scalar_type nsym_norm_error_1=0.0;
+ scalar_type nsym_norm_error_max=0.0;
+ if (tlist[iterm].var1==tlist[iterm].var2){
+ model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
+ gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
+ nsym_norm_error_euc
+ = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
+ nsym_norm_error_1
+ = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
+ nsym_norm_error_max
+ = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
+ }
-
- //print matrix if the size is small
- if(rvc1[0].size()<8){
- std::cout << "RHS Stiffness Matrix: \n";
- std::cout << "------------------------\n";
- for(size_type i=0; i < rvc1[iterm].size(); ++i){
- std::cout << "[";
- for(size_type j=0; j < rvc1[iterm].size(); ++j){
- std::cout << fdSM(i,j) << " ";
- }
- std::cout << "]\n";
- }
- std::cout << "Analytical Stiffness Matrix: \n";
- std::cout << "------------------------\n";
- for(size_type i=0; i < rvc1[iterm].size(); ++i){
- std::cout << "[";
- for(size_type j=0; j < rvc1[iterm].size(); ++j){
- std::cout << matl[iterm](i,j) << " ";
- }
- std::cout << "]\n";
- }
- std::cout << "Vector U: \n";
- std::cout << "------------------------\n";
- for(size_type i=0; i < rvc1[iterm].size(); ++i){
- std::cout << "[";
- std::cout << md.real_variable(tlist[iterm].var2)[i] << " ";
- std::cout << "]\n";
- }
- }
- std::cout
- << "\n\nfinite diff test error_norm_eucl: " << norm_error_euc << "%\n"
- << "finite diff test error_norm1: " << norm_error_1 << "%\n"
- << "finite diff test error_max_norm: " << norm_error_max << "%\n\n\n";
-
- if (tlist[iterm].var1==tlist[iterm].var2){
- std::cout
- << "Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<< "%\n"
- << "Nonsymmetrical test error_norm1: " << nsym_norm_error_1 << "%\n"
- << "Nonsymmetrical test error_max_norm: " << nsym_norm_error_max << "%"
- << std::endl;
- }
- }
++=======
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
+
+ void virtual_brick::check_stiffness_matrix_and_rhs
+ (const model &md, size_type s,
+ const model::termlist& tlist,
+ const model::varnamelist &vl,
+ const model::varnamelist &dl,
+ const model::mimlist &mims,
+ model::real_matlist &matl,
+ model::real_veclist &rvc1,
+ model::real_veclist &rvc2,
+ size_type rg,
+ const scalar_type TINY) const {
+ std::cout<<"******Verifying stiffnesses of *******"<<std::endl;
+ std::cout<<"*** "<<brick_name()<<std::endl;
+
+ //Build the index for the corresponding RHS
+ std::map<std::string,size_type> rhs_index;
+ for(size_type iterm=0;iterm<matl.size();iterm++)
+ if (tlist[iterm].var1==tlist[iterm].var2) rhs_index[tlist[iterm].var1]=iterm;
+
+ if (rhs_index.size()==0){
+ GMM_WARNING0("*** cannot verify stiffness for this brick***");
+ return;
}
++<<<<<<< HEAD
+>>>>>>> upstream
+
++=======
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2,
+ rg, model::BUILD_MATRIX);
+ for(size_type iterm=0;iterm<matl.size();iterm++){
+
+ std::cout<<std::endl;
+ std::cout<<" Stiffness["<<tlist[iterm].var1
+ <<","<<tlist[iterm].var2<<"]:"<<std::endl;
+ if (md.real_variable(tlist[iterm].var1).size()==0)
+ {
+ std::cout<<" "<<tlist[iterm].var1<<" has zero size. Skipping this term"<<std::endl;
+ continue;
+ }
+ if (md.real_variable(tlist[iterm].var2).size()==0)
+ {
+ std::cout<<" "<<tlist[iterm].var2<<" has zero size. Skipping this term"<<std::endl;
+ continue;
+ }
+
+ model_real_sparse_matrix SM(matl[iterm]);
+ gmm::fill(rvc1[rhs_index[tlist[iterm].var1]], 0.0);
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2,
+ rg, model::BUILD_RHS);
+ if (gmm::mat_euclidean_norm(matl[iterm])<1e-12){
+ std::cout<<" The assembled matrix is nearly zero, skipping."<<std::endl;
+ continue;
+ }
+ model_real_plain_vector RHS0(rvc1[rhs_index[tlist[iterm].var1]]);
+
+ //finite difference stiffness
+ model_real_sparse_matrix fdSM(matl[iterm].nrows(),matl[iterm].ncols());
+ model_real_plain_vector&U = md.set_real_variable(tlist[iterm].var2);
+ model_real_plain_vector& RHS1 =rvc1[rhs_index[tlist[iterm].var1]];
+ for (size_type j=0; j < matl[iterm].ncols(); j++){
+ U[j]+=TINY;
+ gmm::fill(RHS1, 0.0);
+ asm_real_tangent_terms(md, s, vl, dl, mims, matl, rvc1, rvc2,
+ rg, model::BUILD_RHS);
+ for (size_type i=0;i<matl[iterm].nrows();i++)
+ fdSM(i,j) = (RHS0[i]-RHS1[i])/TINY;
+ U[j]-=TINY;
+ }
+
+ model_real_sparse_matrix diffSM(matl[iterm].nrows(),matl[iterm].ncols());
+ gmm::add(SM,gmm::scaled(fdSM,-1.0),diffSM);
+ scalar_type norm_error_euc
+ = gmm::mat_euclidean_norm(diffSM)/gmm::mat_euclidean_norm(SM)*100;
+ scalar_type norm_error_1
+ = gmm::mat_norm1(diffSM)/gmm::mat_norm1(SM)*100;
+ scalar_type norm_error_max
+ = gmm::mat_maxnorm(diffSM)/gmm::mat_maxnorm(SM)*100;
+
+ //checking symmetry of diagonal terms
+ scalar_type nsym_norm_error_euc=0.0;
+ scalar_type nsym_norm_error_1=0.0;
+ scalar_type nsym_norm_error_max=0.0;
+ if (tlist[iterm].var1==tlist[iterm].var2){
+ model_real_sparse_matrix diffSMtransposed(matl[iterm].nrows(),matl[iterm].ncols());
+ gmm::add(gmm::transposed(fdSM),gmm::scaled(fdSM,-1.0),diffSMtransposed);
+ nsym_norm_error_euc
+ = gmm::mat_euclidean_norm(diffSMtransposed)/gmm::mat_euclidean_norm(fdSM)*100;
+ nsym_norm_error_1
+ = gmm::mat_norm1(diffSMtransposed)/gmm::mat_norm1(fdSM)*100;
+ nsym_norm_error_max
+ = gmm::mat_maxnorm(diffSMtransposed)/gmm::mat_maxnorm(fdSM)*100;
+ }
+
+ //print matrix if the size is small
+ if(rvc1[0].size()<8){
+ std::cout << "RHS Stiffness Matrix: \n";
+ std::cout << "------------------------\n";
+ for(size_type i=0; i < rvc1[iterm].size(); ++i){
+ std::cout << "[";
+ for(size_type j=0; j < rvc1[iterm].size(); ++j){
+ std::cout << fdSM(i,j) << " ";
+ }
+ std::cout << "]\n";
+ }
+ std::cout << "Analytical Stiffness Matrix: \n";
+ std::cout << "------------------------\n";
+ for(size_type i=0; i < rvc1[iterm].size(); ++i){
+ std::cout << "[";
+ for(size_type j=0; j < rvc1[iterm].size(); ++j){
+ std::cout << matl[iterm](i,j) << " ";
+ }
+ std::cout << "]\n";
+ }
+ std::cout << "Vector U: \n";
+ std::cout << "------------------------\n";
+ for(size_type i=0; i < rvc1[iterm].size(); ++i){
+ std::cout << "[";
+ std::cout << md.real_variable(tlist[iterm].var2)[i] << " ";
+ std::cout << "]\n";
+ }
+ }
+ std::cout
+ << "\n\nfinite diff test error_norm_eucl: " << norm_error_euc << "%\n"
+ << "finite diff test error_norm1: " << norm_error_1 << "%\n"
+ << "finite diff test error_max_norm: " << norm_error_max << "%\n\n\n";
+
+ if (tlist[iterm].var1==tlist[iterm].var2){
+ std::cout
+ << "Nonsymmetrical test error_norm_eucl: "<< nsym_norm_error_euc<< "%\n"
+ << "Nonsymmetrical test error_norm1: " << nsym_norm_error_1 << "%\n"
+ << "Nonsymmetrical test error_max_norm: " << nsym_norm_error_max << "%"
+ << std::endl;
+ }
+ }
+ }
+
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
// ----------------------------------------------------------------------
//
// Generic elliptic brick
@@@ -2753,22 -2512,7 +2949,26 @@@
asm_normal_source_term(vecl[0], mim, mf_u, *mf_data, A, rg);
else
asm_homogeneous_normal_source_term(vecl[0], mim, mf_u, A, rg);
+<<<<<<< HEAD
++<<<<<<< HEAD
+
+ }
+=======
+>>>>>>> upstream
+
+ virtual scalar_type asm_real_pseudo_potential(const model &md, size_type,
+ const model::varnamelist &vl,
+ const model::varnamelist &,
+ const model::mimlist &,
+ model::real_matlist &,
+ model::real_veclist &vecl,
+ model::real_veclist &,
+ size_type) const {
+ const model_real_plain_vector &u = md.real_variable(vl[0]);
+ return -gmm::vect_sp(vecl[0], u);
++=======
+ md.add_external_load(ib, gmm::vect_norm1(vecl[0]));
++>>>>>>> upstream/4.2.1_beta1_svn4453_dfsg
}
virtual scalar_type asm_real_pseudo_potential(const model &md, size_type,
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/getfem.git
More information about the debian-science-commits
mailing list