[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Thu, 6 Sep 2018 14:13:47 -0400 (EDT) |
branch: code-cleanup
commit 703515faa16b6b2325840444ed0e26f9f628608a
Author: Konstantinos Poulios <address@hidden>
Date: Thu Sep 6 20:13:37 2018 +0200
Typo fixes, white space and constness changes
---
interface/src/gf_asm.cc | 39 +-
src/bgeot_convex_ref.cc | 13 +-
src/bgeot_node_tab.cc | 2 +-
src/getfem/bgeot_convex_ref.h | 2 +-
src/getfem/bgeot_geometric_trans.h | 2 -
src/getfem/getfem_generic_assembly_tree.h | 4 +-
src/getfem/getfem_mesh.h | 2 +-
src/getfem/getfem_models.h | 25 +-
src/getfem_fem.cc | 30 +-
src/getfem_generic_assembly_compile_and_exec.cc | 708 ++++++++++++------------
src/getfem_generic_assembly_semantic.cc | 4 +-
src/getfem_generic_assembly_workspace.cc | 4 +-
src/getfem_global_function.cc | 3 +-
src/getfem_mesh.cc | 6 +-
src/getfem_models.cc | 23 +-
src/getfem_projected_fem.cc | 63 ++-
16 files changed, 473 insertions(+), 457 deletions(-)
diff --git a/interface/src/gf_asm.cc b/interface/src/gf_asm.cc
index 5792ecf..20d97d9 100644
--- a/interface/src/gf_asm.cc
+++ b/interface/src/gf_asm.cc
@@ -743,7 +743,7 @@ void gf_asm(getfemint::mexargs_in& m_in,
getfemint::mexargs_out& m_out) {
Performs the generic assembly of `expression` with the integration
method `mim` on the mesh region of index `region` (-1 means all
- the element of the mesh). The same mesh should be shared by
+ elements of the mesh). The same mesh should be shared by
the integration method and all the finite element methods or
mesh_im_data corresponding to the variables.
@@ -752,28 +752,29 @@ void gf_asm(getfemint::mexargs_in& m_in,
getfemint::mexargs_out& m_out) {
tangent (matrix) (order = 2) is to be computed.
`model` is an optional parameter allowing to take into account
- all variables and data of a model. Optionnally, for the integration
- on the product of two domains, a secondary domain of the model can
- be specified after a 'Secondary_domain' string.
-
- The variables and constant (data) are listed after the
- region number (or optionally the model).
- For each variable/constant, first the variable/constant
- name should be given (as it is referred in the assembly string), then
- 1 if it is a variable or 0 for a constant, then the finite element
- method if it is a fem variable/constant or the mesh_im_data if it is
- data defined on integration points, and the vector representing
- the value of the variable/constant. It is possible to give an arbitrary
- number of variable/constant. The difference between a variable and a
- constant is that automatic differentiation is done with respect to
- variables only (see GetFEM++ user documentation). Test functions are
- only available for variables, not for constants.
+ all variables and data of a model. Note that all enabled variables
+ of the model will occupy space in the returned vector/matrix
+ corresponding to their degrees of freedom in the global system, even
+ if they are not present in `expression`.
+
+ The variables and constants (data) are listed after the region number
+ (or optionally the model).
+ For each variable/constant, a name must be given first (as it is
+ referred in the assembly string), then an integer equal to 1 or 0
+ is expected respectively for declaring a variable or a constant,
+ then the finite element method if it is a fem variable/constant or
+ the mesh_im_data if it is data defined on integration points, and
+ the vector representing the value of the variable/constant.
+ It is possible to give an arbitrary number of variable/constant.
+ The difference between a variable and a constant is that test
+ functions are only available for variables, not for constants.
Note that if several variables are given, the assembly of the
tangent matrix/residual vector will be done considering the order
in the call of the function (the degrees of freedom of the first
- variable, then of the second, and so on). If a model is provided,
- all degrees of freedom of the model will be counted first.
+ variable, then of the second one, and so on). If a model is provided,
+ all degrees of freedom of the model will be counted first, even if
+ some of the model variables do not appear in `expression`.
For example, the L2 norm of a vector field "u" can be computed with::
diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index cc3dc59..e1cf5b4 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -113,7 +113,7 @@ namespace bgeot {
size_type simplexified_tab(pconvex_structure cvs, size_type **tab);
static void simplexify_convex(bgeot::convex_of_reference *cvr,
- mesh_structure &m) {
+ mesh_structure &m) {
pconvex_structure cvs = cvr->structure();
m.clear();
auto basic_cvs = basic_structure(cvs);
@@ -196,10 +196,6 @@ namespace bgeot {
// .push_back(psimplexified_convex);
}
- bool convex_of_reference::is_basic() const {
- return auto_basic;
- }
-
/* should be called on the basic_convex_ref */
const mesh_structure* convex_of_reference::simplexified_convex() const {
GMM_ASSERT1(auto_basic,
@@ -247,6 +243,7 @@ namespace bgeot {
for (; it != ite; e += *it, ++it) r = std::max(r, -(*it));
return std::max(r, e);
}
+
scalar_type is_in_face(short_type f, const base_node &pt) const {
// return zero if pt is in the face of the convex
// negative if the point is on the side of the face where the element is
@@ -413,10 +410,11 @@ namespace bgeot {
else
return gmm::vect_sp(normals_[f], pt) - sqrt(2.)/2.;
}
+
scalar_type is_in(const base_node& pt) const {
// return a negative number if pt is in the convex
scalar_type r = is_in_face(0, pt);
- for (short_type i = 1; i < 5; ++i) r = std::max(r, is_in_face(i, pt));
+ for (short_type f = 1; f < 5; ++f) r = std::max(r, is_in_face(f, pt));
return r;
}
@@ -693,6 +691,7 @@ namespace bgeot {
class equilateral_simplex_of_ref_ : public convex_of_reference {
public:
scalar_type is_in(const base_node &pt) const {
+ GMM_ASSERT1(pt.size() == cvs->dim(), "Dimension does not match");
scalar_type d(0);
for (size_type f = 0; f < normals().size(); ++f) {
const base_node &x0 = (f ? convex<base_node>::points()[f-1]
@@ -702,7 +701,9 @@ namespace bgeot {
}
return d;
}
+
scalar_type is_in_face(short_type f, const base_node &pt) const {
+ GMM_ASSERT1(pt.size() == cvs->dim(), "Dimension does not match");
const base_node &x0 = (f ? convex<base_node>::points()[f-1]
: convex<base_node>::points().back());
return gmm::vect_sp(pt-x0, normals()[f]);
diff --git a/src/bgeot_node_tab.cc b/src/bgeot_node_tab.cc
index 580dce5..6f85a72 100644
--- a/src/bgeot_node_tab.cc
+++ b/src/bgeot_node_tab.cc
@@ -85,7 +85,7 @@ namespace bgeot {
GMM_ASSERT1(false, "Problem in node structure !!");
}
- void node_tab::clear(void) {
+ void node_tab::clear() {
dal::dynamic_tas<base_node>::clear();
sorters = std::vector<sorter>();
max_radius = scalar_type(1e-60);
diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h
index b3ec6ef..7eef34f 100644
--- a/src/getfem/bgeot_convex_ref.h
+++ b/src/getfem/bgeot_convex_ref.h
@@ -108,7 +108,7 @@ namespace bgeot {
* positive in the other side.
*/
virtual scalar_type is_in_face(short_type, const base_node &) const =0;
- bool is_basic() const;
+ bool is_basic() const { return auto_basic; }
/// return the normal vector for each face.
const std::vector<base_small_vector> &normals() const
{ return normals_; }
diff --git a/src/getfem/bgeot_geometric_trans.h
b/src/getfem/bgeot_geometric_trans.h
index a002a2e..cc4a7a1 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -159,8 +159,6 @@ namespace bgeot {
template<class CONT> base_node transform(const base_node &pt,
const CONT &PTAB) const;
base_node transform(const base_node &pt, const base_matrix &G) const;
- /** Compute the gradient at point x, pc is resized to [nb_points() x dim()]
- if the transformation is linear, x is not used at all */
size_type complexity() const { return complexity_; }
virtual ~geometric_trans()
{ DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Geometric transformation"); }
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index 9f54f21..a2c9413 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -454,7 +454,9 @@ namespace getfem {
ga_tree &operator =(const ga_tree &tree) {
clear(); secondary_domain = tree.secondary_domain;
- if (tree.root) copy_node(tree.root,nullptr,root); return *this;
+ if (tree.root)
+ copy_node(tree.root,nullptr,root);
+ return *this;
}
~ga_tree() { clear(); }
diff --git a/src/getfem/getfem_mesh.h b/src/getfem/getfem_mesh.h
index 0c86cfa..1748596 100644
--- a/src/getfem/getfem_mesh.h
+++ b/src/getfem/getfem_mesh.h
@@ -645,7 +645,7 @@ namespace getfem {
*/
mesh_region APIDECL
inner_faces_of_mesh(const mesh &m,
- mesh_region mr = mesh_region::all_convexes());
+ 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/getfem_models.h b/src/getfem/getfem_models.h
index 11b397e..1b71bcc 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -1561,19 +1561,19 @@ namespace getfem {
size_type APIDECL add_linear_term
(model &md, const mesh_im &mim, const std::string &expr,
size_type region = size_type(-1), bool is_sym = false,
- bool is_coercive = false, std::string brickname = "",
+ bool is_coercive = false, const std::string &brickname = "",
bool return_if_nonlin = false);
inline size_type APIDECL add_linear_generic_assembly_brick
(model &md, const mesh_im &mim, const std::string &expr,
size_type region = size_type(-1), bool is_sym = false,
- bool is_coercive = false, std::string brickname = "",
+ bool is_coercive = false, const std::string &brickname = "",
bool return_if_nonlin = false) {
return add_linear_term(md, mim, expr, region, is_sym,
is_coercive, brickname, return_if_nonlin);
}
- /** Add a nonlinear term given by the weak form language expression `expr`
+ /** Add a nonlinear term given by the weak form language expression `expr`
which will be assembled in region `region` and with the integration
method `mim`.
The expression can describe a potential or a weak form. Second order
@@ -1593,7 +1593,7 @@ namespace getfem {
inline size_type APIDECL add_nonlinear_generic_assembly_brick
(model &md, const mesh_im &mim, const std::string &expr,
size_type region = size_type(-1), bool is_sym = false,
- bool is_coercive = false, std::string brickname = "") {
+ bool is_coercive = false, const std::string &brickname = "") {
return add_nonlinear_term(md, mim, expr, region,
is_sym, is_coercive, brickname);
}
@@ -1609,14 +1609,16 @@ namespace getfem {
*/
size_type APIDECL add_source_term
(model &md, const mesh_im &mim, const std::string &expr,
- size_type region = size_type(-1), std::string brickname = "",
- std::string directvarname = std::string(),
+ size_type region = size_type(-1),
+ const std::string &brickname = std::string(),
+ const std::string &directvarname = std::string(),
const std::string &directdataname = std::string(),
bool return_if_nonlin = false);
inline size_type APIDECL add_source_term_generic_assembly_brick
(model &md, const mesh_im &mim, const std::string &expr,
- size_type region = size_type(-1), std::string brickname = "",
- std::string directvarname = std::string(),
+ size_type region = size_type(-1),
+ const std::string &brickname = std::string(),
+ const std::string &directvarname = std::string(),
const std::string &directdataname = std::string(),
bool return_if_nonlin = false) {
return add_source_term(md, mim, expr, region, brickname,
@@ -1632,8 +1634,8 @@ namespace getfem {
size_type APIDECL add_linear_twodomain_term
(model &md, const mesh_im &mim, const std::string &expr,
size_type region, const std::string &secondary_domain,
- bool is_sym = false, bool is_coercive = false, std::string brickname = "",
- bool return_if_nonlin = false);
+ bool is_sym = false, bool is_coercive = false,
+ const std::string &brickname = "", bool return_if_nonlin = false);
/** Adds a nonlinear term given by a weak form language expression like
``add_nonlinear_term`` function but for an integration on a direct
@@ -1656,7 +1658,8 @@ namespace getfem {
size_type APIDECL add_twodomain_source_term
(model &md, const mesh_im &mim, const std::string &expr,
size_type region, const std::string &secondary_domain,
- std::string brickname = "", std::string directvarname = std::string(),
+ const std::string &brickname = std::string(),
+ const std::string &directvarname = std::string(),
const std::string &directdataname = std::string(),
bool return_if_nonlin = false);
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 3b45e9e..263da5c 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1605,7 +1605,7 @@ namespace getfem {
}
/* ******************************************************************** */
- /* P1 element with a bubble base fonction on a face */
+ /* P1 element with a bubble base function on a face */
/* ******************************************************************** */
struct P1_wabbfoaf_ : public PK_fem_ {
@@ -1640,7 +1640,7 @@ namespace getfem {
}
/* ******************************************************************** */
- /* Element RT0 on the simplexes.
*/
+ /* Element RT0 on the simplexes. */
/* ******************************************************************** */
struct P1_RT0_ : public fem<base_poly> {
@@ -1680,7 +1680,7 @@ namespace getfem {
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) M(i, i) *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
- GMM_WARNING2("RT0 : The normal orientation may be not correct");
+ GMM_WARNING2("RT0 : The normal orientation may be incorrect");
}
}
@@ -1737,7 +1737,7 @@ namespace getfem {
/* ******************************************************************** */
- /* Element RT0 on parallelepideds.
*/
+ /* Element RT0 on parallelepideds. */
/* ******************************************************************** */
struct P1_RT0Q_ : public fem<base_poly> {
@@ -1777,7 +1777,7 @@ namespace getfem {
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) M(i, i) *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
- GMM_WARNING2("RT0Q : The normal orientation may be not correct");
+ GMM_WARNING2("RT0Q : The normal orientation may be incorrect");
}
}
@@ -1834,7 +1834,7 @@ namespace getfem {
/* ******************************************************************** */
- /* Nedelec Element.
*/
+ /* Nedelec Element. */
/* ******************************************************************** */
struct P1_nedelec_ : public fem<base_poly> {
@@ -1949,7 +1949,7 @@ namespace getfem {
/* ******************************************************************** */
- /* P1 element with a bubble base fonction on a face : type lagrange
*/
+ /* P1 element with a bubble base function on a face : type lagrange */
/* ******************************************************************** */
struct P1_wabbfoafla_ : public PK_fem_
@@ -1981,7 +1981,7 @@ namespace getfem {
/* ******************************************************************** */
- /* PK Gauss-Lobatto element on the segment
*/
+ /* PK Gauss-Lobatto element on the segment */
/* ******************************************************************** */
static const double fem_coef_gausslob_1[4] =
@@ -3227,7 +3227,7 @@ namespace getfem {
}
/* ******************************************************************** */
- /* Hermite element on the triangle
*/
+ /* Hermite element on the triangle */
/* ******************************************************************** */
struct hermite_triangle__ : public fem<base_poly> {
@@ -3300,7 +3300,7 @@ namespace getfem {
}
/* ******************************************************************** */
- /* Hermite element on the tetrahedron
*/
+ /* Hermite element on the tetrahedron */
/* ******************************************************************** */
struct hermite_tetrahedron__ : public fem<base_poly> {
@@ -3471,7 +3471,7 @@ namespace getfem {
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) n *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
- GMM_WARNING2("Argyris : The normal orientation may be not correct");
+ GMM_WARNING2("Argyris : The normal orientation may be incorrect");
gmm::mult(K, n, v);
const bgeot::base_tensor &t = pfp->grad(i);
for (unsigned j = 0; j < 21; ++j)
@@ -3615,7 +3615,7 @@ namespace getfem {
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) n *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
- GMM_WARNING2("Morley : The normal orientation may be not correct");
+ GMM_WARNING2("Morley : The normal orientation may be incorrect");
gmm::mult(K, n, v);
const bgeot::base_tensor &t = pfp->grad(i);
for (unsigned j = 0; j < 6; ++j)
@@ -3673,7 +3673,7 @@ namespace getfem {
}
/* ******************************************************************** */
- /* DISCONTINUOUS PK
*/
+ /* DISCONTINUOUS PK */
/* ******************************************************************** */
struct PK_discont_ : public PK_fem_ {
@@ -3724,7 +3724,7 @@ namespace getfem {
/* ******************************************************************** */
- /* PK element with a bubble base fonction
*/
+ /* PK element with a bubble base function */
/* ******************************************************************** */
struct PK_with_cubic_bubble_ : public PK_fem_ {
@@ -3771,7 +3771,7 @@ namespace getfem {
}
/* ******************************************************************** */
- /* classical fem
*/
+ /* classical fem */
/* ******************************************************************** */
static pfem classical_fem_(bgeot::pgeometric_trans pgt,
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 3b02b25..1cd22df 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -780,7 +780,7 @@ namespace getfem {
virtual int exec() {
GA_DEBUG_INFO("Instruction: value of test functions");
if (qdim == 1) {
- GA_DEBUG_ASSERT(t.size() == Z.size(), "Wrong size for base vector");
+ GA_DEBUG_ASSERT(t.size() == Z.size(), "Wrong size for base vector");
std::copy(Z.begin(), Z.end(), t.begin());
} else {
size_type target_dim = Z.sizes()[1];
@@ -3796,7 +3796,7 @@ namespace getfem {
if (inin.pt_type) {
if (cv != size_type(-1)) {
inin.m->points_of_convex(cv, inin.G);
- inin.ctx.change((inin.m)->trans_of_convex(cv),
+ inin.ctx.change((inin.m)->trans_of_convex(cv),
0, P_ref, inin.G, cv, face_num);
inin.has_ctx = true;
if (face_num != short_type(-1)) {
@@ -3888,7 +3888,7 @@ namespace getfem {
gic.init(m.points_of_convex(adj_face.cv), gpc.pgt2);
size_type first_ind = pai->ind_first_point_on_face(f);
const bgeot::stored_point_tab
- &spt = *(pai->pintegration_points());
+ &spt = *(pai->pintegration_points());
base_matrix G;
m.points_of_convex(cv, G);
fem_interpolation_context ctx_x(gpc.pgt1, 0, spt[0], G, cv, f);
@@ -3898,8 +3898,8 @@ namespace getfem {
ctx_x.set_xref(spt[first_ind+i]);
bool converged = true;
gic.invert(ctx_x.xreal(), P_ref[i], converged);
- bool is_in = (gpc.pgt2->convex_ref()->is_in(P_ref[i]) < 1E-4);
- GMM_ASSERT1(is_in && converged,"Geometric transformation "
+ bool is_in = (gpc.pgt2->convex_ref()->is_in(P_ref[i]) < 1E-4);
+ GMM_ASSERT1(is_in && converged,"Geometric transformation "
"inversion has failed in neighbour
transformation");
}
pspt = store_point_tab(P_ref);
@@ -3953,7 +3953,7 @@ namespace getfem {
}
}
GA_DEBUG_INFO("Instruction: end of call neighbour interpolate "
- "transformation");
+ "transformation");
return 0;
}
ga_instruction_neighbour_transformation_call
@@ -4772,17 +4772,17 @@ namespace getfem {
pctx1 = &(gis.ctx);
const std::string &intn1 = pnode->interpolate_name_test1;
if (intn1.size()) {
- if (workspace.secondary_domain_exists(intn1)) {
- pctx1 = &(rmi.secondary_domain_infos.ctx);
- } else {
- tensor_to_adapt = true;
- pctx1 = &(rmi.interpolate_infos[intn1].ctx);
- if (workspace.variable_group_exists(pnode->name_test1)) {
- ga_instruction_set::variable_group_info &vgi =
- rmi.interpolate_infos[intn1].groups_info[pnode->name_test1];
- mfg1 = &(vgi.mf);
- mf1 = 0;
- }
+ if (workspace.secondary_domain_exists(intn1)) {
+ pctx1 = &(rmi.secondary_domain_infos.ctx);
+ } else {
+ tensor_to_adapt = true;
+ pctx1 = &(rmi.interpolate_infos[intn1].ctx);
+ if (workspace.variable_group_exists(pnode->name_test1)) {
+ ga_instruction_set::variable_group_info &vgi =
+ rmi.interpolate_infos[intn1].groups_info[pnode->name_test1];
+ mfg1 = &(vgi.mf);
+ mf1 = 0;
+ }
}
}
}
@@ -4792,18 +4792,18 @@ namespace getfem {
pctx2 = &(gis.ctx);
const std::string &intn2 = pnode->interpolate_name_test2;
if (intn2.size()) {
- if (workspace.secondary_domain_exists(intn2)) {
- pctx2 = &(rmi.secondary_domain_infos.ctx);
- } else {
- tensor_to_adapt = true;
- pctx2 = &(rmi.interpolate_infos[intn2].ctx);
- if (workspace.variable_group_exists(pnode->name_test2)) {
- ga_instruction_set::variable_group_info &vgi =
- rmi.interpolate_infos[intn2].groups_info[pnode->name_test2];
- mfg2 = &(vgi.mf);
- mf2 = 0;
- }
- }
+ if (workspace.secondary_domain_exists(intn2)) {
+ pctx2 = &(rmi.secondary_domain_infos.ctx);
+ } else {
+ tensor_to_adapt = true;
+ pctx2 = &(rmi.interpolate_infos[intn2].ctx);
+ if (workspace.variable_group_exists(pnode->name_test2)) {
+ ga_instruction_set::variable_group_info &vgi =
+ rmi.interpolate_infos[intn2].groups_info[pnode->name_test2];
+ mfg2 = &(vgi.mf);
+ mf2 = 0;
+ }
+ }
}
}
}
@@ -4850,17 +4850,15 @@ namespace getfem {
// Optimization: detects if an equivalent node has already been compiled
pnode->t.set_to_original();
if (rmi.node_list.find(pnode->hash_value) != rmi.node_list.end()) {
- std::list<pga_tree_node> &node_list = rmi.node_list[pnode->hash_value];
- for (std::list<pga_tree_node>::iterator it = node_list.begin();
- it != node_list.end(); ++it) {
- // cout << "found potential equivalent nodes ";
- // ga_print_node(pnode, cout);
- // cout << " and "; ga_print_node(*it, cout); cout << endl;
- if (sub_tree_are_equal(pnode, *it, workspace, 1)) {
- pnode->t.set_to_copy((*it)->t);
+ for (pga_tree_node &pnode1 : rmi.node_list[pnode->hash_value]) {
+ // cout << "found potential equivalent nodes ";
+ // ga_print_node(pnode, cout);
+ // cout << " and "; ga_print_node(pnode1, cout); cout << endl;
+ if (sub_tree_are_equal(pnode, pnode1, workspace, 1)) {
+ pnode->t.set_to_copy(pnode1->t);
return;
}
- if (sub_tree_are_equal(pnode, *it, workspace, 2)) {
+ if (sub_tree_are_equal(pnode, pnode1, workspace, 2)) {
// cout << "confirmed with transpose" << endl;
if (pnode->nb_test_functions() == 2) {
if (pgai) { // resize instruction if needed
@@ -4869,17 +4867,17 @@ namespace getfem {
else { rmi.instructions.push_back(std::move(pgai)); }
}
pgai = std::make_shared<ga_instruction_transpose_test>
- (pnode->tensor(), (*it)->tensor());
+ (pnode->tensor(), pnode1->tensor());
rmi.instructions.push_back(std::move(pgai));
} else {
- pnode->t.set_to_copy((*it)->t);
+ pnode->t.set_to_copy(pnode1->t);
}
return;
}
std::stringstream ss;
ss << "Detected wrong equivalent nodes: ";
ga_print_node(pnode, ss);
- ss << " and "; ga_print_node(*it, ss);
+ ss << " and "; ga_print_node(pnode1, ss);
ss << " (no problem, but hash values could be adapted) " << endl;
GMM_TRACE2(ss.str());
}
@@ -5023,19 +5021,19 @@ namespace getfem {
case GA_NODE_SECONDARY_DOMAIN_X:
case GA_NODE_SECONDARY_DOMAIN_NORMAL:
{
- GMM_ASSERT1(!function_case,
- "No use of Secondary_domain is allowed in functions");
- auto psd = workspace.secondary_domain(pnode->interpolate_name);
- size_type sddim = psd->mim().linked_mesh().dim();
- if (pnode->tensor().size() != sddim)
- pnode->init_vector_tensor(sddim);
- if (pnode->node_type == GA_NODE_SECONDARY_DOMAIN_X)
- pgai = std::make_shared<ga_instruction_X>
- (pnode->tensor(), rmi.secondary_domain_infos.ctx);
- else if (pnode->node_type == GA_NODE_SECONDARY_DOMAIN_NORMAL)
- pgai = std::make_shared<ga_instruction_copy_Normal>
- (pnode->tensor(), rmi.secondary_domain_infos.Normal);
- rmi.instructions.push_back(std::move(pgai));
+ GMM_ASSERT1(!function_case,
+ "No use of Secondary_domain is allowed in functions");
+ auto psd = workspace.secondary_domain(pnode->interpolate_name);
+ size_type sddim = psd->mim().linked_mesh().dim();
+ if (pnode->tensor().size() != sddim)
+ pnode->init_vector_tensor(sddim);
+ if (pnode->node_type == GA_NODE_SECONDARY_DOMAIN_X)
+ pgai = std::make_shared<ga_instruction_X>
+ (pnode->tensor(), rmi.secondary_domain_infos.ctx);
+ else if (pnode->node_type == GA_NODE_SECONDARY_DOMAIN_NORMAL)
+ pgai = std::make_shared<ga_instruction_copy_Normal>
+ (pnode->tensor(), rmi.secondary_domain_infos.Normal);
+ rmi.instructions.push_back(std::move(pgai));
}
break;
@@ -5319,16 +5317,16 @@ namespace getfem {
case GA_NODE_SECONDARY_DOMAIN_VAL: case GA_NODE_SECONDARY_DOMAIN_GRAD:
case GA_NODE_SECONDARY_DOMAIN_HESS: case GA_NODE_SECONDARY_DOMAIN_DIVERG:
{
- GMM_ASSERT1(!function_case, "internal error");
+ GMM_ASSERT1(!function_case, "internal error");
const mesh_fem *mf = workspace.associated_mf(pnode->name);
const im_data *imd = workspace.associated_im_data(pnode->name);
- const std::string &intn = pnode->interpolate_name;
- auto &sdi = rmi.secondary_domain_infos;
+ const std::string &intn = pnode->interpolate_name;
+ auto &sdi = rmi.secondary_domain_infos;
+
+ fem_interpolation_context *pctx = &(sdi.ctx);
+ papprox_integration pai = sdi.pai;
+ psecondary_domain psd = workspace.secondary_domain(intn);
- fem_interpolation_context *pctx = &(sdi.ctx);
- papprox_integration pai = sdi.pai;
- psecondary_domain psd = workspace.secondary_domain(intn);
-
if (imd) {
pgai = std::make_shared<ga_instruction_extract_local_im_data>
(pnode->tensor(), *imd, workspace.value(pnode->name),
@@ -5340,7 +5338,7 @@ namespace getfem {
"The finite element of variable " << pnode->name <<
" has to be defined on the same mesh than the "
"integration method or interpolation used on the "
- "secondary domain");
+ "secondary domain");
// An instruction for extracting local dofs of the variable.
if (sdi.local_dofs.count(pnode->name) == 0) {
@@ -5379,7 +5377,7 @@ namespace getfem {
}
break;
case GA_NODE_SECONDARY_DOMAIN_GRAD:
- case GA_NODE_SECONDARY_DOMAIN_DIVERG:
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG:
if (sdi.grad.count(mf) == 0 ||
!(if_hierarchy.is_compatible(rmi.grad_hierarchy[mf]))) {
rmi.grad_hierarchy[mf].push_back(if_hierarchy);
@@ -5742,14 +5740,14 @@ namespace getfem {
case GA_NODE_SECONDARY_DOMAIN_HESS_TEST:
case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
{
- GMM_ASSERT1(!function_case, "internal error");
- const mesh_fem *mf = workspace.associated_mf(pnode->name);
- const std::string &intn = pnode->interpolate_name;
- auto &sdi = rmi.secondary_domain_infos;
-
- fem_interpolation_context *pctx = &(sdi.ctx);
- papprox_integration pai = sdi.pai;
- psecondary_domain psd = workspace.secondary_domain(intn);
+ GMM_ASSERT1(!function_case, "internal error");
+ const mesh_fem *mf = workspace.associated_mf(pnode->name);
+ const std::string &intn = pnode->interpolate_name;
+ auto &sdi = rmi.secondary_domain_infos;
+
+ fem_interpolation_context *pctx = &(sdi.ctx);
+ papprox_integration pai = sdi.pai;
+ psecondary_domain psd = workspace.secondary_domain(intn);
if (mf) {
GMM_ASSERT1(&(mf->linked_mesh()) == &(psd->mim().linked_mesh()),
"The finite element of variable " << pnode->name <<
@@ -5779,7 +5777,7 @@ namespace getfem {
}
break;
case GA_NODE_SECONDARY_DOMAIN_GRAD_TEST:
- case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
+ case GA_NODE_SECONDARY_DOMAIN_DIVERG_TEST:
if (sdi.grad.count(mf) == 0 ||
!(if_hierarchy.is_compatible(rmi.grad_hierarchy[mf]))) {
rmi.grad_hierarchy[mf].push_back(if_hierarchy);
@@ -6699,9 +6697,9 @@ namespace getfem {
// Compile tree
// cout << "Will compile "; ga_print_node(root, cout); cout <<
endl;
- psecondary_domain psd(0);
- if (added_tree->secondary_domain.size())
- psd = workspace.secondary_domain(added_tree->secondary_domain);
+ psecondary_domain psd(0);
+ if (added_tree->secondary_domain.size())
+ psd = workspace.secondary_domain(added_tree->secondary_domain);
ga_instruction_set::region_mim rm(td.mim, td.rg, psd);
ga_instruction_set::region_mim_instructions &rmi
= gis.whole_instructions[rm];
@@ -6745,8 +6743,8 @@ namespace getfem {
if (mf) {
const std::string &intn1 = root->interpolate_name_test1;
const gmm::sub_interval *Ir = 0, *In = 0;
- bool secondary = intn1.size() &&
- workspace.secondary_domain_exists(intn1);
+ bool secondary = intn1.size() &&
+ workspace.secondary_domain_exists(intn1);
if (intn1.size() && !secondary &&
workspace.variable_group_exists(root->name_test1)) {
ga_instruction_set::variable_group_info &vgi =
@@ -6762,11 +6760,11 @@ namespace getfem {
}
fem_interpolation_context &ctx
= intn1.size() ?
- (secondary ? rmi.secondary_domain_infos.ctx
- : rmi.interpolate_infos[intn1].ctx) : gis.ctx;
+ (secondary ? rmi.secondary_domain_infos.ctx
+ : rmi.interpolate_infos[intn1].ctx) : gis.ctx;
bool interpolate
= (!intn1.empty() && intn1.compare("neighbour_elt")!=0 &&
- !secondary);
+ !secondary);
pgai = std::make_shared<ga_instruction_fem_vector_assembly>
(root->tensor(), workspace.unreduced_vector(),
workspace.assembled_vector(), ctx, *Ir, *In, mf, mfg,
@@ -6790,22 +6788,22 @@ namespace getfem {
const std::string &intn1 = root->interpolate_name_test1;
const std::string &intn2 = root->interpolate_name_test2;
bool secondary1 = intn1.size() &&
- workspace.secondary_domain_exists(intn1);
+ workspace.secondary_domain_exists(intn1);
bool secondary2 = intn2.size() &&
- workspace.secondary_domain_exists(intn2);
- fem_interpolation_context &ctx1
+ workspace.secondary_domain_exists(intn2);
+ fem_interpolation_context &ctx1
= intn1.size() ?
- (secondary1 ? rmi.secondary_domain_infos.ctx
- : rmi.interpolate_infos[intn1].ctx) : gis.ctx;
- fem_interpolation_context &ctx2
+ (secondary1 ? rmi.secondary_domain_infos.ctx
+ : rmi.interpolate_infos[intn1].ctx) : gis.ctx;
+ fem_interpolation_context &ctx2
= intn2.size() ?
- (secondary2 ? rmi.secondary_domain_infos.ctx
- : rmi.interpolate_infos[intn2].ctx) : gis.ctx;
- bool interpolate
- = (!intn1.empty() && intn1.compare("neighbour_elt")!=0 &&
- !secondary1) ||
- (!intn2.empty() && intn2.compare("neighbour_elt")!=0 &&
- !secondary2);
+ (secondary2 ? rmi.secondary_domain_infos.ctx
+ : rmi.interpolate_infos[intn2].ctx) : gis.ctx;
+ bool interpolate
+ = (!intn1.empty() && intn1.compare("neighbour_elt")!=0 &&
+ !secondary1) ||
+ (!intn2.empty() && intn2.compare("neighbour_elt")!=0 &&
+ !secondary2);
add_interval_to_gis(workspace, root->name_test1, gis);
add_interval_to_gis(workspace, root->name_test2, gis);
@@ -7053,272 +7051,272 @@ namespace getfem {
if (!psd) { // standard integration on a single domain
- const mesh_region ®ion = *(instr.first.region());
-
- // iteration on elements (or faces of elements)
- size_type old_cv = size_type(-1);
- bgeot::pgeometric_trans pgt = 0, pgt_old = 0;
- pintegration_method pim = 0;
- papprox_integration pai = 0;
- bgeot::pstored_point_tab pspt = 0, old_pspt = 0;
- bgeot::pgeotrans_precomp pgp = 0;
- bool first_gp = true;
- for (getfem::mr_visitor v(region, m, true); !v.finished(); ++v) {
- if (mim.convex_index().is_in(v.cv())) {
- // cout << "proceed with elt " << v.cv() << " face " << v.f()<<endl;
- if (v.cv() != old_cv) {
- pgt = m.trans_of_convex(v.cv());
- pim = mim.int_method_of_element(v.cv());
- m.points_of_convex(v.cv(), G1);
-
- if (pim->type() == IM_NONE) continue;
- GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods "
- "cannot be used in high level generic assembly");
- pai = pim->approx_method();
- pspt = pai->pintegration_points();
- if (pspt->size()) {
- if (pgp && gis.pai == pai && pgt_old == pgt) {
- gis.ctx.change(pgp, 0, 0, G1, v.cv(), v.f());
- } else {
- if (pai->is_built_on_the_fly()) {
- gis.ctx.change(pgt, 0, (*pspt)[0], G1, v.cv(), v.f());
- pgp = 0;
- } else {
- pgp = gis.gp_pool(pgt, pspt);
- gis.ctx.change(pgp, 0, 0, G1, v.cv(), v.f());
- }
- pgt_old = pgt; gis.pai = pai;
- }
- if (gis.need_elt_size)
- gis.elt_size = convex_radius_estimate(pgt, G1)*scalar_type(2);
- }
- old_cv = v.cv();
- } else {
- if (pim->type() == IM_NONE) continue;
- gis.ctx.set_face_num(v.f());
- }
- if (pspt != old_pspt) { first_gp = true; old_pspt = pspt; }
- if (pspt->size()) {
- // iterations on Gauss points
- size_type first_ind = 0;
- if (v.f() != short_type(-1)) {
- gis.nbpt = pai->nb_points_on_face(v.f());
- first_ind = pai->ind_first_point_on_face(v.f());
- } else {
- gis.nbpt = pai->nb_points_on_convex();
- }
- for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
- if (pgp) gis.ctx.set_ii(first_ind+gis.ipt);
- else gis.ctx.set_xref((*pspt)[first_ind+gis.ipt]);
- if (gis.ipt == 0 || !(pgt->is_linear())) {
- J1 = gis.ctx.J();
- // Computation of unit normal vector in case of a boundary
- if (v.f() != short_type(-1)) {
- gis.Normal.resize(G1.nrows());
- un.resize(pgt->dim());
- gmm::copy(pgt->normals()[v.f()], un);
- gmm::mult(gis.ctx.B(), un, gis.Normal);
- scalar_type nup = gmm::vect_norm2(gis.Normal);
- J1 *= nup;
- gmm::scale(gis.Normal, 1.0/nup);
- gmm::clean(gis.Normal, 1e-13);
- } else gis.Normal.resize(0);
- }
- auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
- gis.coeff = J1 * ipt_coeff;
- bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
- workspace.include_empty_int_points());
- if (!enable_ipt) gis.coeff = scalar_type(0);
- if (first_gp) {
- for (size_type j=0; j < gilb.size(); ++j) j+=gilb[j]->exec();
- first_gp = false;
- }
- if (gis.ipt == 0) {
- for (size_type j=0; j < gile.size(); ++j) j+=gile[j]->exec();
- }
- if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1) {
- for (size_type j=0; j < gil.size(); ++j) j+=gil[j]->exec();
- }
- GA_DEBUG_INFO("");
- }
- }
- }
- }
- GA_DEBUG_INFO("-----------------------------");
-
+ const mesh_region ®ion = *(instr.first.region());
+
+ // iteration on elements (or faces of elements)
+ size_type old_cv = size_type(-1);
+ bgeot::pgeometric_trans pgt = 0, pgt_old = 0;
+ pintegration_method pim = 0;
+ papprox_integration pai = 0;
+ bgeot::pstored_point_tab pspt = 0, old_pspt = 0;
+ bgeot::pgeotrans_precomp pgp = 0;
+ bool first_gp = true;
+ for (getfem::mr_visitor v(region, m, true); !v.finished(); ++v) {
+ if (mim.convex_index().is_in(v.cv())) {
+ // cout << "proceed with elt " << v.cv() << " face " <<
v.f()<<endl;
+ if (v.cv() != old_cv) {
+ pgt = m.trans_of_convex(v.cv());
+ pim = mim.int_method_of_element(v.cv());
+ m.points_of_convex(v.cv(), G1);
+
+ if (pim->type() == IM_NONE) continue;
+ GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods "
+ "cannot be used in high level generic assembly");
+ pai = pim->approx_method();
+ pspt = pai->pintegration_points();
+ if (pspt->size()) {
+ if (pgp && gis.pai == pai && pgt_old == pgt) {
+ gis.ctx.change(pgp, 0, 0, G1, v.cv(), v.f());
+ } else {
+ if (pai->is_built_on_the_fly()) {
+ gis.ctx.change(pgt, 0, (*pspt)[0], G1, v.cv(), v.f());
+ pgp = 0;
+ } else {
+ pgp = gis.gp_pool(pgt, pspt);
+ gis.ctx.change(pgp, 0, 0, G1, v.cv(), v.f());
+ }
+ pgt_old = pgt; gis.pai = pai;
+ }
+ if (gis.need_elt_size)
+ gis.elt_size = convex_radius_estimate(pgt,
G1)*scalar_type(2);
+ }
+ old_cv = v.cv();
+ } else {
+ if (pim->type() == IM_NONE) continue;
+ gis.ctx.set_face_num(v.f());
+ }
+ if (pspt != old_pspt) { first_gp = true; old_pspt = pspt; }
+ if (pspt->size()) {
+ // iterations on Gauss points
+ size_type first_ind = 0;
+ if (v.f() != short_type(-1)) {
+ gis.nbpt = pai->nb_points_on_face(v.f());
+ first_ind = pai->ind_first_point_on_face(v.f());
+ } else {
+ gis.nbpt = pai->nb_points_on_convex();
+ }
+ for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
+ if (pgp) gis.ctx.set_ii(first_ind+gis.ipt);
+ else gis.ctx.set_xref((*pspt)[first_ind+gis.ipt]);
+ if (gis.ipt == 0 || !(pgt->is_linear())) {
+ J1 = gis.ctx.J();
+ // Computation of unit normal vector in case of a boundary
+ if (v.f() != short_type(-1)) {
+ gis.Normal.resize(G1.nrows());
+ un.resize(pgt->dim());
+ gmm::copy(pgt->normals()[v.f()], un);
+ gmm::mult(gis.ctx.B(), un, gis.Normal);
+ scalar_type nup = gmm::vect_norm2(gis.Normal);
+ J1 *= nup;
+ gmm::scale(gis.Normal, 1.0/nup);
+ gmm::clean(gis.Normal, 1e-13);
+ } else gis.Normal.resize(0);
+ }
+ auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
+ gis.coeff = J1 * ipt_coeff;
+ bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
+ workspace.include_empty_int_points());
+ if (!enable_ipt) gis.coeff = scalar_type(0);
+ if (first_gp) {
+ for (size_type j=0; j < gilb.size(); ++j) j+=gilb[j]->exec();
+ first_gp = false;
+ }
+ if (gis.ipt == 0) {
+ for (size_type j=0; j < gile.size(); ++j) j+=gile[j]->exec();
+ }
+ if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1) {
+ for (size_type j=0; j < gil.size(); ++j) j+=gil[j]->exec();
+ }
+ GA_DEBUG_INFO("");
+ }
+ }
+ }
+ }
+ GA_DEBUG_INFO("-----------------------------");
+
} else { // Integration on the product of two domains (secondary domain)
- auto &sdi = instr.second.secondary_domain_infos;
- const mesh_region ®ion1 = *(instr.first.region());
-
- // iteration on elements (or faces of elements)
- size_type old_cv1=size_type(-1), old_cv2=size_type(-1);
- size_type nbpt1 = 0, nbpt2 = 0;
- bgeot::pgeometric_trans pgt1 = 0, pgt1_old = 0, pgt2 = 0, pgt2_old = 0;
- pintegration_method pim1 = 0, pim2 = 0;
- papprox_integration pai1 = 0, pai2 = 0;
- bgeot::pstored_point_tab pspt1=0, old_pspt1=0, pspt2=0, old_pspt2=0;
- bgeot::pgeotrans_precomp pgp1 = 0, pgp2 = 0;
- bool first_gp = true;
- for (getfem::mr_visitor v1(region1, m, true); !v1.finished(); ++v1) {
- if (mim.convex_index().is_in(v1.cv())) {
- // cout << "proceed with elt " << v1.cv()<<" face " << v1.f()<<endl;
- if (v1.cv() != old_cv1) {
- pgt1 = m.trans_of_convex(v1.cv());
- pim1 = mim.int_method_of_element(v1.cv());
- m.points_of_convex(v1.cv(), G1);
-
- if (pim1->type() == IM_NONE) continue;
- GMM_ASSERT1(pim1->type() == IM_APPROX, "Sorry, exact methods "
- "cannot be used in high level generic assembly");
- pai1 = pim1->approx_method();
- pspt1 = pai1->pintegration_points();
- if (pspt1->size()) {
- if (pgp1 && gis.pai == pai1 && pgt1_old == pgt1) {
- gis.ctx.change(pgp1, 0, 0, G1, v1.cv(), v1.f());
- } else {
- if (pai1->is_built_on_the_fly()) {
- gis.ctx.change(pgt1, 0, (*pspt1)[0], G1, v1.cv(), v1.f());
- pgp1 = 0;
- } else {
- pgp1 = gis.gp_pool(pgt1, pspt1);
- gis.ctx.change(pgp1, 0, 0, G1, v1.cv(), v1.f());
- }
- pgt1_old = pgt1; gis.pai = pai1;
- }
- if (gis.need_elt_size)
- gis.elt_size = convex_radius_estimate(pgt1,G1)*scalar_type(2);
- }
- old_cv1 = v1.cv();
- } else {
- if (pim1->type() == IM_NONE) continue;
- gis.ctx.set_face_num(v1.f());
- }
- if (pspt1 != old_pspt1) { first_gp = true; old_pspt1 = pspt1; }
- if (pspt1->size()) {
- // iterations on Gauss points
- size_type first_ind1 = 0;
- if (v1.f() != short_type(-1)) {
- nbpt1 = pai1->nb_points_on_face(v1.f());
- first_ind1 = pai1->ind_first_point_on_face(v1.f());
- } else {
- nbpt1 = pai1->nb_points_on_convex();
- }
-
- const mesh &m2 = psd->mim().linked_mesh();
- const mesh_region ®ion2 = psd->give_region(m, v1.cv(), v1.f());
- for (getfem::mr_visitor v2(region2, m2, true);
- !v2.finished(); ++v2) {
- if (v2.cv() != old_cv2) {
- pgt2 = m2.trans_of_convex(v2.cv());
- pim2 = psd->mim().int_method_of_element(v2.cv());
- m2.points_of_convex(v2.cv(), G2);
-
- if (pim2->type() == IM_NONE) continue;
- GMM_ASSERT1(pim2->type() == IM_APPROX, "Sorry, exact methods "
- "cannot be used in high level generic assembly");
- pai2 = pim2->approx_method();
- pspt2 = pai2->pintegration_points();
- if (pspt2->size()) {
- if (pgp2 && sdi.pai == pai2 && pgt2_old == pgt2) {
- sdi.ctx.change(pgp2, 0, 0, G2, v2.cv(), v2.f());
- } else {
- if (pai2->is_built_on_the_fly()) {
- sdi.ctx.change(pgt2, 0, (*pspt2)[0], G2,v2.cv(),v2.f());
- pgp2 = 0;
- } else {
- pgp2 = gis.gp_pool(pgt2, pspt2);
- sdi.ctx.change(pgp2, 0, 0, G2, v2.cv(), v2.f());
- }
- pgt2_old = pgt2; sdi.pai = pai2;
- }
- }
- old_cv2 = v2.cv();
- } else {
- if (pim2->type() == IM_NONE) continue;
- sdi.ctx.set_face_num(v2.f());
- }
- if (pspt2 != old_pspt2) { first_gp = true; old_pspt2 = pspt2; }
- if (pspt2->size()) {
- // iterations on Gauss points
- size_type first_ind2 = 0;
- if (v2.f() != short_type(-1)) {
- nbpt2 = pai2->nb_points_on_face(v2.f());
- first_ind2 = pai2->ind_first_point_on_face(v2.f());
- } else {
- nbpt2 = gis.nbpt = pai2->nb_points_on_convex();
- }
- gis.nbpt = nbpt1 * nbpt2;
- gis.ipt = 0;
- for (size_type ipt1=0; ipt1 < nbpt1; ++ipt1) {
- for (size_type ipt2=0; ipt2 < nbpt2; ++ipt2, ++(gis.ipt)) {
-
- if (pgp1) gis.ctx.set_ii(first_ind1+ipt1);
- else gis.ctx.set_xref((*pspt1)[first_ind1+ipt1]);
- if (pgp2) sdi.ctx.set_ii(first_ind2+ipt2);
- else sdi.ctx.set_xref((*pspt2)[first_ind2+ipt2]);
-
- if (gis.ipt == 0 || !(pgt1->is_linear())) {
- J1 = gis.ctx.J();
- if (v1.f() != short_type(-1)) {
- gis.Normal.resize(G1.nrows());
- un.resize(pgt1->dim());
- gmm::copy(pgt1->normals()[v1.f()], un);
- gmm::mult(gis.ctx.B(), un, gis.Normal);
- scalar_type nup = gmm::vect_norm2(gis.Normal);
- J1 *= nup;
- gmm::scale(gis.Normal, 1.0/nup);
- gmm::clean(gis.Normal, 1e-13);
- } else gis.Normal.resize(0);
- }
-
- if (gis.ipt == 0 || !(pgt2->is_linear())) {
- J2 = sdi.ctx.J();
- if (v2.f() != short_type(-1)) {
- sdi.Normal.resize(G2.nrows());
- un.resize(pgt2->dim());
- gmm::copy(pgt2->normals()[v2.f()], un);
- gmm::mult(sdi.ctx.B(), un, sdi.Normal);
- scalar_type nup = gmm::vect_norm2(sdi.Normal);
- J2 *= nup;
- gmm::scale(sdi.Normal, 1.0/nup);
- gmm::clean(sdi.Normal, 1e-13);
- } else sdi.Normal.resize(0);
- }
-
- auto ipt_coeff = pai1->coeff(first_ind1+ipt1)
- * pai2->coeff(first_ind2+ipt2);
- gis.coeff = J1 * J2 * ipt_coeff;
- bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
- workspace.include_empty_int_points());
- if (!enable_ipt) gis.coeff = scalar_type(0);
-
- if (first_gp) {
- for (size_type j=0; j < gilb.size(); ++j)
- j+=gilb[j]->exec();
- first_gp = false;
- }
- if (gis.ipt == 0) {
- for (size_type j=0; j < gile.size(); ++j)
- j+=gile[j]->exec();
- }
- if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1) {
- for (size_type j=0; j < gil.size(); ++j)
- j+=gil[j]->exec();
- }
- GA_DEBUG_INFO("");
- }
- }
- }
- }
- }
- }
- }
- GA_DEBUG_INFO("-----------------------------");
+ auto &sdi = instr.second.secondary_domain_infos;
+ const mesh_region ®ion1 = *(instr.first.region());
+
+ // iteration on elements (or faces of elements)
+ size_type old_cv1=size_type(-1), old_cv2=size_type(-1);
+ size_type nbpt1 = 0, nbpt2 = 0;
+ bgeot::pgeometric_trans pgt1 = 0, pgt1_old = 0, pgt2 = 0, pgt2_old = 0;
+ pintegration_method pim1 = 0, pim2 = 0;
+ papprox_integration pai1 = 0, pai2 = 0;
+ bgeot::pstored_point_tab pspt1=0, old_pspt1=0, pspt2=0, old_pspt2=0;
+ bgeot::pgeotrans_precomp pgp1 = 0, pgp2 = 0;
+ bool first_gp = true;
+ for (getfem::mr_visitor v1(region1, m, true); !v1.finished(); ++v1) {
+ if (mim.convex_index().is_in(v1.cv())) {
+ // cout << "proceed with elt " << v1.cv()<<" face " <<
v1.f()<<endl;
+ if (v1.cv() != old_cv1) {
+ pgt1 = m.trans_of_convex(v1.cv());
+ pim1 = mim.int_method_of_element(v1.cv());
+ m.points_of_convex(v1.cv(), G1);
+
+ if (pim1->type() == IM_NONE) continue;
+ GMM_ASSERT1(pim1->type() == IM_APPROX, "Sorry, exact methods "
+ "cannot be used in high level generic assembly");
+ pai1 = pim1->approx_method();
+ pspt1 = pai1->pintegration_points();
+ if (pspt1->size()) {
+ if (pgp1 && gis.pai == pai1 && pgt1_old == pgt1) {
+ gis.ctx.change(pgp1, 0, 0, G1, v1.cv(), v1.f());
+ } else {
+ if (pai1->is_built_on_the_fly()) {
+ gis.ctx.change(pgt1, 0, (*pspt1)[0], G1, v1.cv(), v1.f());
+ pgp1 = 0;
+ } else {
+ pgp1 = gis.gp_pool(pgt1, pspt1);
+ gis.ctx.change(pgp1, 0, 0, G1, v1.cv(), v1.f());
+ }
+ pgt1_old = pgt1; gis.pai = pai1;
+ }
+ if (gis.need_elt_size)
+ gis.elt_size =
convex_radius_estimate(pgt1,G1)*scalar_type(2);
+ }
+ old_cv1 = v1.cv();
+ } else {
+ if (pim1->type() == IM_NONE) continue;
+ gis.ctx.set_face_num(v1.f());
+ }
+ if (pspt1 != old_pspt1) { first_gp = true; old_pspt1 = pspt1; }
+ if (pspt1->size()) {
+ // iterations on Gauss points
+ size_type first_ind1 = 0;
+ if (v1.f() != short_type(-1)) {
+ nbpt1 = pai1->nb_points_on_face(v1.f());
+ first_ind1 = pai1->ind_first_point_on_face(v1.f());
+ } else {
+ nbpt1 = pai1->nb_points_on_convex();
+ }
+
+ const mesh &m2 = psd->mim().linked_mesh();
+ const mesh_region ®ion2 = psd->give_region(m, v1.cv(),
v1.f());
+ for (getfem::mr_visitor v2(region2, m2, true);
+ !v2.finished(); ++v2) {
+ if (v2.cv() != old_cv2) {
+ pgt2 = m2.trans_of_convex(v2.cv());
+ pim2 = psd->mim().int_method_of_element(v2.cv());
+ m2.points_of_convex(v2.cv(), G2);
+
+ if (pim2->type() == IM_NONE) continue;
+ GMM_ASSERT1(pim2->type() == IM_APPROX, "Sorry, exact methods
"
+ "cannot be used in high level generic assembly");
+ pai2 = pim2->approx_method();
+ pspt2 = pai2->pintegration_points();
+ if (pspt2->size()) {
+ if (pgp2 && sdi.pai == pai2 && pgt2_old == pgt2) {
+ sdi.ctx.change(pgp2, 0, 0, G2, v2.cv(), v2.f());
+ } else {
+ if (pai2->is_built_on_the_fly()) {
+ sdi.ctx.change(pgt2, 0, (*pspt2)[0],
G2,v2.cv(),v2.f());
+ pgp2 = 0;
+ } else {
+ pgp2 = gis.gp_pool(pgt2, pspt2);
+ sdi.ctx.change(pgp2, 0, 0, G2, v2.cv(), v2.f());
+ }
+ pgt2_old = pgt2; sdi.pai = pai2;
+ }
+ }
+ old_cv2 = v2.cv();
+ } else {
+ if (pim2->type() == IM_NONE) continue;
+ sdi.ctx.set_face_num(v2.f());
+ }
+ if (pspt2 != old_pspt2) { first_gp = true; old_pspt2 = pspt2; }
+ if (pspt2->size()) {
+ // iterations on Gauss points
+ size_type first_ind2 = 0;
+ if (v2.f() != short_type(-1)) {
+ nbpt2 = pai2->nb_points_on_face(v2.f());
+ first_ind2 = pai2->ind_first_point_on_face(v2.f());
+ } else {
+ nbpt2 = gis.nbpt = pai2->nb_points_on_convex();
+ }
+ gis.nbpt = nbpt1 * nbpt2;
+ gis.ipt = 0;
+ for (size_type ipt1=0; ipt1 < nbpt1; ++ipt1) {
+ for (size_type ipt2=0; ipt2 < nbpt2; ++ipt2, ++(gis.ipt)) {
+
+ if (pgp1) gis.ctx.set_ii(first_ind1+ipt1);
+ else gis.ctx.set_xref((*pspt1)[first_ind1+ipt1]);
+ if (pgp2) sdi.ctx.set_ii(first_ind2+ipt2);
+ else sdi.ctx.set_xref((*pspt2)[first_ind2+ipt2]);
+
+ if (gis.ipt == 0 || !(pgt1->is_linear())) {
+ J1 = gis.ctx.J();
+ if (v1.f() != short_type(-1)) {
+ gis.Normal.resize(G1.nrows());
+ un.resize(pgt1->dim());
+ gmm::copy(pgt1->normals()[v1.f()], un);
+ gmm::mult(gis.ctx.B(), un, gis.Normal);
+ scalar_type nup = gmm::vect_norm2(gis.Normal);
+ J1 *= nup;
+ gmm::scale(gis.Normal, 1.0/nup);
+ gmm::clean(gis.Normal, 1e-13);
+ } else gis.Normal.resize(0);
+ }
+
+ if (gis.ipt == 0 || !(pgt2->is_linear())) {
+ J2 = sdi.ctx.J();
+ if (v2.f() != short_type(-1)) {
+ sdi.Normal.resize(G2.nrows());
+ un.resize(pgt2->dim());
+ gmm::copy(pgt2->normals()[v2.f()], un);
+ gmm::mult(sdi.ctx.B(), un, sdi.Normal);
+ scalar_type nup = gmm::vect_norm2(sdi.Normal);
+ J2 *= nup;
+ gmm::scale(sdi.Normal, 1.0/nup);
+ gmm::clean(sdi.Normal, 1e-13);
+ } else sdi.Normal.resize(0);
+ }
+
+ auto ipt_coeff = pai1->coeff(first_ind1+ipt1)
+ * pai2->coeff(first_ind2+ipt2);
+ gis.coeff = J1 * J2 * ipt_coeff;
+ bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
+ workspace.include_empty_int_points());
+ if (!enable_ipt) gis.coeff = scalar_type(0);
+
+ if (first_gp) {
+ for (size_type j=0; j < gilb.size(); ++j)
+ j+=gilb[j]->exec();
+ first_gp = false;
+ }
+ if (gis.ipt == 0) {
+ for (size_type j=0; j < gile.size(); ++j)
+ j+=gile[j]->exec();
+ }
+ if (enable_ipt || gis.ipt == 0 || gis.ipt == gis.nbpt-1)
{
+ for (size_type j=0; j < gil.size(); ++j)
+ j+=gil[j]->exec();
+ }
+ GA_DEBUG_INFO("");
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ GA_DEBUG_INFO("-----------------------------");
}
-
+
}
-
+
for (const std::string &t : gis.transformations)
workspace.interpolate_transformation(t)->finalize();
}
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 41a3850..54be6d9 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -2575,7 +2575,7 @@ namespace getfem {
size_type ref_elt_dim,
bool eval_fixed_size,
bool ignore_X, int option) {
- // cout << "Begin semantic anaylsis" << endl;
+ // cout << "Begin semantic analysis" << endl;
GMM_ASSERT1(predef_operators_nonlinear_elasticity_initialized &&
predef_operators_plasticity_initialized &&
predef_operators_contact_initialized, "Internal error");
@@ -2596,7 +2596,7 @@ namespace getfem {
tree.clear();
}
ga_valid_operand(tree.root);
- // cout << "end of semantic anaylsis" << endl;
+ // cout << "end of semantic analysis" << endl;
}
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index e750ab8..af872ec 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -316,7 +316,7 @@ namespace getfem {
pinterpolate_transformation
ga_workspace::interpolate_transformation(const std::string &name) const {
- auto it = transformations.find(name);
+ auto it = transformations.find(name);
if (it != transformations.end()) return it->second;
if (md && md->interpolate_transformation_exists(name))
return md->interpolate_transformation(name);
@@ -336,7 +336,7 @@ namespace getfem {
pelementary_transformation
ga_workspace::elementary_transformation(const std::string &name) const {
- auto it = elem_transformations.find(name);
+ auto it = elem_transformations.find(name);
if (it != elem_transformations.end()) return it->second;
if (md && md->elementary_transformation_exists(name))
return md->elementary_transformation(name);
diff --git a/src/getfem_global_function.cc b/src/getfem_global_function.cc
index c0fca22..1e140d0 100644
--- a/src/getfem_global_function.cc
+++ b/src/getfem_global_function.cc
@@ -507,7 +507,8 @@ namespace getfem {
return res;
}
- base_matrix crack_singular_xy_function::hess(scalar_type x, scalar_type y)
const {
+ base_matrix
+ crack_singular_xy_function::hess(scalar_type x, scalar_type y) const {
scalar_type sgny = (y < 0 ? -1.0 : 1.0);
scalar_type r = sqrt(x*x + y*y);
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index 5c6ebd8..5a79eed 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -815,7 +815,7 @@ namespace getfem {
between the two neighbour elements. Try to minimize the number of
elements.
*/
- mesh_region inner_faces_of_mesh(const mesh &m, mesh_region mr) {
+ mesh_region inner_faces_of_mesh(const mesh &m, const mesh_region &mr) {
mesh_region mrr;
mr.from_mesh(m);
mr.error_if_not_convexes();
@@ -926,8 +926,8 @@ namespace getfem {
out.clear();
size_type nbpt = in.points().index().last()+1;
GMM_ASSERT1(nbpt == in.points().index().card(),
- "please optimize the mesh before using "
- "it as a base for prismatic mesh");
+ "please call the optimize_structure() method before using "
+ "the mesh as a base for prismatic mesh");
size_type nb_layers_total = nb_layers * degree;
for (const base_node &inpt : in.points()) {
std::copy(inpt.begin(), inpt.end(), pt.begin());
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index afa4f71..294d173 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -3363,7 +3363,7 @@ model_complex_plain_vector &
};
static bool check_compatibility_vl_test(model &md,
- const model::varnamelist vl_test) {
+ const model::varnamelist vl_test) {
model::varnamelist org;
for (size_type i = 0; i < vl_test.size(); ++i) {
if (md.is_affine_dependent_variable(vl_test[i]))
@@ -3377,7 +3377,7 @@ model_complex_plain_vector &
size_type add_source_term_
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
- std::string brickname, std::string directvarname,
+ const std::string &brickname, std::string directvarname,
const std::string &directdataname, bool return_if_nonlin,
const std::string &secondary_domain) {
@@ -3414,19 +3414,19 @@ model_complex_plain_vector &
size_type add_source_term
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
- std::string brickname, std::string directvarname,
+ const std::string &brickname, const std::string &directvarname,
const std::string &directdataname, bool return_if_nonlin) {
return add_source_term_(md, mim, expr, region, brickname, directvarname,
- directdataname, return_if_nonlin, "");
+ directdataname, return_if_nonlin, "");
}
size_type add_twodomain_source_term
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
const std::string &secondary_domain,
- std::string brickname, std::string directvarname,
+ const std::string &brickname, const std::string &directvarname,
const std::string &directdataname, bool return_if_nonlin) {
return add_source_term_(md, mim, expr, region, brickname, directvarname,
- directdataname, return_if_nonlin, secondary_domain);
+ directdataname, return_if_nonlin,
secondary_domain);
}
// ----------------------------------------------------------------------
@@ -3538,7 +3538,7 @@ model_complex_plain_vector &
size_type add_linear_term_
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
- bool is_sym, bool is_coercive, std::string brickname,
+ bool is_sym, bool is_coercive, const std::string &brickname,
bool return_if_nonlin, const std::string &secondary_domain) {
ga_workspace workspace(md, true);
@@ -3580,7 +3580,7 @@ model_complex_plain_vector &
size_type add_linear_term
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
- bool is_sym, bool is_coercive, std::string brickname,
+ bool is_sym, bool is_coercive, const std::string &brickname,
bool return_if_nonlin) {
return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
brickname, return_if_nonlin, "");
@@ -3589,7 +3589,7 @@ model_complex_plain_vector &
size_type add_linear_twodomain_term
(model &md, const mesh_im &mim, const std::string &expr, size_type region,
const std::string &secondary_domain, bool is_sym, bool is_coercive,
- std::string brickname, bool return_if_nonlin) {
+ const std::string &brickname, bool return_if_nonlin) {
return add_linear_term_(md, mim, expr, region, is_sym, is_coercive,
brickname, return_if_nonlin, secondary_domain);
}
@@ -5885,6 +5885,11 @@ model_complex_plain_vector &
GMM_ASSERT1(mims.size() == 0, "Explicit matrix need no mesh_im");
GMM_ASSERT1(vl.size() >= 1 && vl.size() <= 2 && dl.size() == 0,
"Wrong number of variables for explicit matrix brick");
+ GMM_ASSERT1(gmm::mat_ncols(rB) == gmm::mat_ncols(matl[0]) &&
+ gmm::mat_nrows(rB) == gmm::mat_nrows(matl[0]),
+ "Explicit matrix brick dimension mismatch ("<<
+ gmm::mat_ncols(rB)<<"x"<<gmm::mat_nrows(rB)<<") != ("<<
+ gmm::mat_ncols(matl[0])<<"x"<<gmm::mat_nrows(matl[0])<<")");
gmm::copy(rB, matl[0]);
}
diff --git a/src/getfem_projected_fem.cc b/src/getfem_projected_fem.cc
index 10d5534..34adf8d 100644
--- a/src/getfem_projected_fem.cc
+++ b/src/getfem_projected_fem.cc
@@ -283,7 +283,7 @@ namespace getfem {
else { // project on convex faces
mesh_region::face_bitset faces = rg_source.faces_of_convex(cv);
if (faces.count() > 0) { // this should rarely be more than one face
- bgeot::vectors_to_base_matrix(G,
mf_source.linked_mesh().points_of_convex(cv));
+ mf_source.linked_mesh().points_of_convex(cv, G);
short_type nbf = mf_source.linked_mesh().nb_faces_of_convex(cv);
for (short_type f = 0; f < nbf; ++f) {
if (faces.test(f)) {
@@ -337,8 +337,8 @@ namespace getfem {
dim_ = dim__;
if (i.is_face()) dim__ = dim_type(dim__ - 1);
GMM_ASSERT1(dim__ < N, "The projection should take place in lower "
- "dimensions than the mesh dimension. Otherwise "
- "use the interpolated_fem object instead.");
+ "dimensions than the mesh dimension. Otherwise "
+ "use the interpolated_fem object instead.");
}
else
GMM_ASSERT1(dim_ == dim__,
@@ -366,7 +366,7 @@ namespace getfem {
if (gppd.iflags) {
// calculate gppd.normal
const bgeot::pgeometric_trans pgt_source =
mf_source.linked_mesh().trans_of_convex(gppd.cv);
- bgeot::vectors_to_base_matrix(G,
mf_source.linked_mesh().points_of_convex(gppd.cv));
+ mf_source.linked_mesh().points_of_convex(gppd.cv, G);
if (gppd.f != short_type(-1))
normal_on_convex_face(pgt_source, G, gppd.f, gppd.ptref,
gppd.normal);
else
@@ -380,15 +380,17 @@ namespace getfem {
if (gppd.f == short_type(-1)) { // convex
size_type nbdof = mf_source.nb_basic_dof_of_element(gppd.cv);
for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) {
- size_type idof =
mf_source.ind_basic_dof_of_element(gppd.cv)[loc_dof];
- if (!(blocked_dofs[idof])) dofs.add(idof);
+ size_type dof =
mf_source.ind_basic_dof_of_element(gppd.cv)[loc_dof];
+ if (!(blocked_dofs[dof]))
+ dofs.add(dof);
}
}
else { // convex face
size_type nbdof =
mf_source.nb_basic_dof_of_face_of_element(gppd.cv, gppd.f);
for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) {
- size_type idof =
mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
- if (!(blocked_dofs[idof])) dofs.add(idof);
+ size_type dof =
mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
+ if (!(blocked_dofs[dof]))
+ dofs.add(dof);
}
}
last_cv = gppd.cv;
@@ -400,17 +402,17 @@ namespace getfem {
e.inddof.resize(dofs.card());
max_dof = std::max(max_dof, dofs.card());
size_type cnt = 0;
- for (dal::bv_visitor idof(dofs); !idof.finished(); ++idof)
- { e.inddof[cnt] = idof; ind_dof[idof] = cnt++; }
+ for (dal::bv_visitor dof(dofs); !dof.finished(); ++dof)
+ { e.inddof[cnt] = dof; ind_dof[dof] = cnt++; }
for (size_type k = 0; k < nb_pts; ++k) {
gausspt_projection_data &gppd = e.gausspt[start_pt + k];
if (gppd.iflags) {
if (gppd.f == short_type(-1)) { // convex
size_type nbdof = mf_source.nb_basic_dof_of_element(gppd.cv);
for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) {
- size_type idof =
mf_source.ind_basic_dof_of_element(gppd.cv)[loc_dof];
- gppd.local_dof[loc_dof] = dofs.is_in(idof) ? ind_dof[idof]
- : size_type(-1);
+ size_type dof =
mf_source.ind_basic_dof_of_element(gppd.cv)[loc_dof];
+ gppd.local_dof[loc_dof] = dofs.is_in(dof) ? ind_dof[dof]
+ : size_type(-1);
}
}
else { // convex face
@@ -420,19 +422,19 @@ namespace getfem {
unsigned rdim = target_dim() / pf->target_dim();
if (rdim == 1)
for (size_type loc_dof = 0; loc_dof < nbdof; ++loc_dof) { //
local dof with respect to the source convex face
- size_type idof =
mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
+ size_type dof =
mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
size_type loc_dof2 = ind_pts_fc[loc_dof]; // local dof with
respect to the source convex
- gppd.local_dof[loc_dof2] = dofs.is_in(idof) ? ind_dof[idof]
- : size_type(-1);
+ gppd.local_dof[loc_dof2] = dofs.is_in(dof) ? ind_dof[dof]
+ : size_type(-1);
}
else
for (size_type ii = 0; ii < nbdof/rdim; ++ii)
for (size_type jj = 0; jj < rdim; ++jj) {
size_type loc_dof = ii*rdim + jj; // local dof with respect
to the source convex face
- size_type idof =
mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
+ size_type dof =
mf_source.ind_basic_dof_of_face_of_element(gppd.cv, gppd.f)[loc_dof];
size_type loc_dof2 = ind_pts_fc[ii]*rdim + jj; // local dof
with respect to the source convex
- gppd.local_dof[loc_dof2] = dofs.is_in(idof) ? ind_dof[idof]
- : size_type(-1);
+ gppd.local_dof[loc_dof2] = dofs.is_in(dof) ? ind_dof[dof]
+ : size_type(-1);
}
}
}
@@ -447,7 +449,7 @@ namespace getfem {
std::fill(dof_types_.begin(), dof_types_.end(),
global_dof(dim()));
- /* ind_dof should be kept full of -1 ( real_base_value and
+ /* ind_dof should be kept filled with -1 ( real_base_value and
grad_base_value expect that)
*/
std::fill(ind_dof.begin(), ind_dof.end(), size_type(-1));
@@ -495,8 +497,7 @@ namespace getfem {
inline void projected_fem::actualize_fictx(pfem pf, size_type cv,
const base_node &ptr) const {
if (fictx_cv != cv) {
- bgeot::vectors_to_base_matrix
- (G, mf_source.linked_mesh().points_of_convex(cv));
+ mf_source.linked_mesh().points_of_convex(cv, G);
fictx = fem_interpolation_context
(mf_source.linked_mesh().trans_of_convex(cv), pf, base_node(), G, cv);
fictx_cv = cv;
@@ -509,7 +510,8 @@ namespace getfem {
std::map<size_type,elt_projection_data>::iterator eit;
eit = elements.find(c.convex_num());
if (eit == elements.end()) {
- mi2[1] = target_dim(); mi2[0] = short_type(0);
+ mi2[1] = target_dim();
+ mi2[0] = short_type(0);
t.adjust_sizes(mi2);
std::fill(t.begin(), t.end(), scalar_type(0));
return;
@@ -517,7 +519,8 @@ namespace getfem {
// GMM_ASSERT1(eit != elements.end(), "Wrong convex number: " <<
c.convex_num());
elt_projection_data &e = eit->second;
- mi2[1] = target_dim(); mi2[0] = short_type(e.nb_dof);
+ mi2[1] = target_dim();
+ mi2[0] = short_type(e.nb_dof);
t.adjust_sizes(mi2);
std::fill(t.begin(), t.end(), scalar_type(0));
if (e.nb_dof == 0) return;
@@ -601,8 +604,10 @@ namespace getfem {
std::map<size_type,elt_projection_data>::iterator eit;
eit = elements.find(c.convex_num());
if (eit == elements.end()) {
- mi3[2] = mf_source.linked_mesh().dim(); mi3[1] = target_dim(); mi3[0] =
short_type(0);
- t.adjust_sizes(mi2);
+ mi3[2] = mf_source.linked_mesh().dim();
+ mi3[1] = target_dim();
+ mi3[0] = short_type(0);
+ t.adjust_sizes(mi3);
std::fill(t.begin(), t.end(), scalar_type(0));
return;
}
@@ -610,7 +615,9 @@ namespace getfem {
elt_projection_data &e = eit->second;
size_type N = mf_source.linked_mesh().dim();
- mi3[2] = short_type(N); mi3[1] = target_dim(); mi3[0] =
short_type(e.nb_dof);
+ mi3[2] = short_type(N);
+ mi3[1] = target_dim();
+ mi3[0] = short_type(e.nb_dof);
t.adjust_sizes(mi3);
std::fill(t.begin(), t.end(), scalar_type(0));
if (e.nb_dof == 0) return;
@@ -735,7 +742,7 @@ namespace getfem {
short_type f;
if (find_a_projected_point(pt, ptref, cv, f)) {
const bgeot::pgeometric_trans pgt =
mf_source.linked_mesh().trans_of_convex(cv);
- bgeot::vectors_to_base_matrix(G,
mf_source.linked_mesh().points_of_convex(cv));
+ mf_source.linked_mesh().points_of_convex(cv, G);
if (f != short_type(-1))
normal_on_convex_face(pgt, G, f, ptref, normal);
else