[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: |
Mon, 11 Jun 2018 05:25:17 -0400 (EDT) |
branch: master
commit 6a5e4f88805ce9fdf31c4ac2fa292b476dfb83b4
Author: Konstantinos Poulios <address@hidden>
Date: Mon Jun 11 11:24:59 2018 +0200
White space and typo fixes
---
.../source/userdoc/model_generic_assembly.rst | 2 +-
interface/tests/matlab/tutorial1.m | 2 +-
src/bgeot_geotrans_inv.cc | 102 ++--
src/getfem/getfem_mesh_im_level_set.h | 2 +-
src/getfem/getfem_models.h | 2 +-
src/getfem_generic_assembly_compile_and_exec.cc | 675 +++++++++++----------
src/getfem_generic_assembly_interpolation.cc | 6 +-
src/getfem_level_set.cc | 30 +-
src/getfem_mesh.cc | 256 ++++----
src/getfem_mesh_im_level_set.cc | 2 -
src/getfem_mesh_region.cc | 254 ++++----
src/getfem_models.cc | 76 ++-
12 files changed, 701 insertions(+), 708 deletions(-)
diff --git a/doc/sphinx/source/userdoc/model_generic_assembly.rst
b/doc/sphinx/source/userdoc/model_generic_assembly.rst
index 43f9cd9..eb355a1 100644
--- a/doc/sphinx/source/userdoc/model_generic_assembly.rst
+++ b/doc/sphinx/source/userdoc/model_generic_assembly.rst
@@ -27,7 +27,7 @@ However, this brick consider that the expression is
nonlinear. This brick is esp
with the same arguments. Conversely, this brick alway assume that the term
corresponding to ``expr`` is linear and the assembly will be performed only
once if the data used do not change. Thus, you have to care that your
expression is indeed linear (affine in fact) with respect to each variable.
Otherwise, the result is of course not guaranted. Source terms in the
expression are taken into account. Still for linear problem, it is possible to
perform the assembly of a sole source term tha [...]
- size_type getfem::add_source(md, mim, expr, region = -1);
+ size_type getfem::add_source_term(md, mim, expr, region = -1);
with again the same arguments except the symmetry and coercivness. This brick
performs the assembly of the corresponding order 1 term (residual vector) and
add it as a right hand side to the problem. The assembly will be performed only
once, so the term should not depend on the variables of the model (but could
depend of course on the constants).
diff --git a/interface/tests/matlab/tutorial1.m
b/interface/tests/matlab/tutorial1.m
index 77db022..3e4ecf2 100644
--- a/interface/tests/matlab/tutorial1.m
+++ b/interface/tests/matlab/tutorial1.m
@@ -26,7 +26,7 @@ mim=gf_mesh_im(m, gf_integ('IM_GAUSS_PARALLELEPIPED(2, 4)'));
% mim=gf_mesh_im(m, gf_integ('IM_EXACT_PARALLELEPIPED(2)')); % not allowed
with the high level generic assembly
border = gf_mesh_get(m,'outer faces');
-gf_mesh_set(m, 'region', 42, border); % create the region (:#42
+gf_mesh_set(m, 'region', 42, border); % create the region (:#42
% the boundary edges appears in red
gf_plot_mesh(m, 'regions', [42], 'vertices','on','convexes','on');
diff --git a/src/bgeot_geotrans_inv.cc b/src/bgeot_geotrans_inv.cc
index ccabaf3..f87e9df 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -28,9 +28,9 @@ using namespace gmm;
using namespace std;
namespace bgeot
-{
+{
bool geotrans_inv_convex::invert(const base_node& n, base_node& n_ref,
- scalar_type IN_EPS) {
+ scalar_type IN_EPS) {
assert(pgt);
n_ref.resize(pgt->structure()->dim());
bool converged = true;
@@ -41,9 +41,9 @@ namespace bgeot
}
}
- bool geotrans_inv_convex::invert(const base_node& n, base_node& n_ref,
- bool &converged,
- scalar_type IN_EPS) {
+ bool geotrans_inv_convex::invert(const base_node& n, base_node& n_ref,
+ bool &converged,
+ scalar_type IN_EPS) {
assert(pgt);
n_ref.resize(pgt->structure()->dim());
converged = true;
@@ -53,24 +53,24 @@ namespace bgeot
}
bool point_in_convex(const geometric_trans &geoTrans,
- const base_node &x,
- scalar_type res,
- scalar_type IN_EPS) {
+ const base_node &x,
+ scalar_type res,
+ scalar_type IN_EPS) {
// Test un peu sev�re peut-�tre en ce qui concerne res.
return (geoTrans.convex_ref()->is_in(x) < IN_EPS) && (res < IN_EPS);
}
/* inversion for linear geometric transformations */
bool geotrans_inv_convex::invert_lin(const base_node& n, base_node& n_ref,
- scalar_type IN_EPS) {
+ scalar_type IN_EPS) {
base_node y(n); for (size_type i=0; i < N; ++i) y[i] -= G(i,0);
mult(transposed(B), y, n_ref);
- y = pgt->transform(n_ref, G);
+ y = pgt->transform(n_ref, G);
add(scaled(n, -1.0), y);
return point_in_convex(*pgt, n_ref, vect_norm2(y), IN_EPS);
}
-
+
void geotrans_inv_convex::update_B() {
if (P != N) {
pgt->compute_K_matrix(G, pc, K);
@@ -85,7 +85,7 @@ namespace bgeot
pgt->compute_K_matrix(G, pc, KT);
gmm::copy(gmm::transposed(KT), K);
gmm::copy(K, B);
- bgeot::lu_inverse(&(*(K.begin())), P); B.swap(K);
+ bgeot::lu_inverse(&(*(K.begin())), P); B.swap(K);
}
}
@@ -97,8 +97,8 @@ namespace bgeot
geotrans_inv_convex &gic;
base_node xreal;
public:
- geotrans_inv_convex_bfgs(geotrans_inv_convex &gic_,
- const base_node &xr) : gic(gic_), xreal(xr) {}
+ geotrans_inv_convex_bfgs(geotrans_inv_convex &gic_,
+ const base_node &xr) : gic(gic_), xreal(xr) {}
scalar_type operator()(const base_node& x) const {
base_node r = gic.pgt->transform(x, gic.G) - xreal;
return gmm::vect_norm2_sqr(r)/2.;
@@ -108,7 +108,7 @@ namespace bgeot
gic.update_B();
base_node r = gic.pgt->transform(x, gic.G) - xreal;
gr.resize(x.size());
- gmm::mult(gmm::transposed(gic.K), r, gr);
+ gmm::mult(gmm::transposed(gic.K), r, gr);
}
};
@@ -122,12 +122,12 @@ namespace bgeot
std::vector<base_node> direct_points_ref;
direct_points.reserve(n_points);
direct_points_ref.reserve(n_points);
-
+
for (auto i : direct_points_indices) {
direct_points.push_back(real_nodes[i]);
direct_points_ref.push_back(reference_nodes[i]);
}
-
+
auto N = direct_points.begin()->size();
auto N_ref = direct_points_ref.begin()->size();
base_matrix K_linear(N, n_points - 1);
@@ -135,13 +135,13 @@ namespace bgeot
K_ref_linear.base_resize(N_ref, n_points - 1);
P_linear = direct_points[0];
P_ref_linear = direct_points_ref[0];
-
+
for (size_type i = 1; i < n_points; ++i) {
add(direct_points[i], scaled(P_linear, -1.0), mat_col(K_linear, i - 1));
add(direct_points_ref[i], scaled(P_ref_linear, -1.0),
- mat_col(K_ref_linear, i - 1));
+ mat_col(K_ref_linear, i - 1));
}
-
+
if (K_linear.nrows() == K_linear.ncols()) {
lu_inverse(K_linear);
copy(transposed(K_linear), B_linear);
@@ -163,79 +163,79 @@ namespace bgeot
}
void project_into_convex(base_node &x, const geometric_trans *pgeo_trans,
- bool project) {
+ bool project) {
if (project == false) return;
-
+
auto dim = pgeo_trans->dim();
-
+
for (auto d = 0; d < dim; ++d) {
if (x[d] < 0.0) x[d] = 0.0;
if (x[d] > 1.0) x[d] = 1.0;
}
auto poriginal_trans = pgeo_trans;
-
+
if (auto ptorus_trans = dynamic_cast<const torus_geom_trans*>(pgeo_trans))
{
poriginal_trans = ptorus_trans->get_original_transformation().get();
}
-
+
auto pbasic_convex_ref = basic_convex_ref(poriginal_trans->convex_ref());
auto nb_simplices = pbasic_convex_ref->simplexified_convex()->nb_convex();
-
+
if (nb_simplices == 1) { // simplex
auto sum_coordinates = 0.0;
-
+
for (auto d = 0; d < dim; ++d) sum_coordinates += x[d];
-
+
if (sum_coordinates > 1.0) scale(x, 1.0 / sum_coordinates);
}
else if ((dim == 3) && (nb_simplices != 4)) { // prism
auto sum_coordinates = x[0] + x[1];
-
+
if (sum_coordinates > 1.0) {
- x[0] /= sum_coordinates;
- x[1] /= sum_coordinates;
+ x[0] /= sum_coordinates;
+ x[1] /= sum_coordinates;
}
}
}
-
+
void find_initial_guess(base_node &x,
- nonlinear_storage_struct &storage,
- const base_node &xreal,
- const base_matrix &G,
- const geometric_trans *pgt,
- scalar_type IN_EPS) {
+ nonlinear_storage_struct &storage,
+ const base_node &xreal,
+ const base_matrix &G,
+ const geometric_trans *pgt,
+ scalar_type IN_EPS) {
storage.x_ref = pgt->geometric_nodes()[0];
-
+
auto res = vect_dist2(mat_col(G, 0), xreal);
double res0;
-
- for (size_type j = 1; j < pgt->nb_points(); ++j) {
+
+ for (size_type j = 1; j < pgt->nb_points(); ++j) {
res0 = vect_dist2(mat_col(G, j), xreal);
if (res > res0) {
- res = res0;
- storage.x_ref = pgt->geometric_nodes()[j];
+ res = res0;
+ storage.x_ref = pgt->geometric_nodes()[j];
}
}
-
+
res0 = std::numeric_limits<scalar_type>::max();
-
+
if (storage.plinearised_structure != nullptr) {
storage.plinearised_structure->invert(xreal, x, IN_EPS);
project_into_convex(x, pgt, storage.project_into_element);
res0 = vect_dist2(pgt->transform(x, G), xreal);
}
-
+
if (res < res0) copy(storage.x_ref, x);
}
-
- /* inversion for non-linear geometric transformations
+
+ /* inversion for non-linear geometric transformations
(Newton on Grad(pgt)(y - pgt(x)) = 0 )
*/
bool geotrans_inv_convex::invert_nonlin(const base_node& xreal,
- base_node& x, scalar_type IN_EPS,
- bool &converged, bool /* throw_except */) {
+ base_node& x, scalar_type IN_EPS,
+ bool &converged, bool /* throw_except */) {
converged = true;
find_initial_guess(x, nonlinear_storage, xreal, G, pgt.get(), IN_EPS);
add(pgt->transform(x, G), scaled(xreal, -1.0), nonlinear_storage.diff);
@@ -254,7 +254,7 @@ namespace bgeot
add(scaled(nonlinear_storage.x_ref, factor), x);
nonlinear_storage.x_real = pgt->transform(x, G);
add(nonlinear_storage.x_real, scaled(xreal, -1.0),
- nonlinear_storage.diff);
+ nonlinear_storage.diff);
factor *= 0.5;
}
else {
@@ -269,7 +269,7 @@ namespace bgeot
project_into_convex(x, pgt.get(),
nonlinear_storage.project_into_element);
nonlinear_storage.x_real = pgt->transform(x, G);
add(nonlinear_storage.x_real, scaled(xreal, -1.0),
- nonlinear_storage.diff);
+ nonlinear_storage.diff);
res = vect_norm2(nonlinear_storage.diff);
++cnt;
}
diff --git a/src/getfem/getfem_mesh_im_level_set.h
b/src/getfem/getfem_mesh_im_level_set.h
index 623d267..0d92830 100644
--- a/src/getfem/getfem_mesh_im_level_set.h
+++ b/src/getfem/getfem_mesh_im_level_set.h
@@ -112,7 +112,7 @@ namespace getfem {
base_singular_pim = sing;
}
- int location(void) const { return integrate_where; }
+ int location() const { return integrate_where; }
size_type memsize() const {
return mesh_im::memsize(); // + ... ;
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 0172c37..a0aeb22 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -2199,7 +2199,7 @@ namespace getfem {
/** Linear elasticity brick ( @f$ \int \sigma(u):\varepsilon(v) @f$ ).
- for isotropic material. Parametrized by the Lam� coefficients
+ for isotropic material. Parametrized by the Lamé coefficients
lambda and mu.
*/
size_type APIDECL add_isotropic_linearized_elasticity_brick
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 27d50db..fd9c330 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -43,7 +43,7 @@ namespace getfem {
return (gpc1.pgt2 < gpc2.pgt2);
return false;
}
-
+
//=========================================================================
// Instructions for compilation: basic optimized operations on tensors
//=========================================================================
@@ -76,7 +76,8 @@ namespace getfem {
papprox_integration &pai_, const fem_interpolation_context &ctx_,
size_type qdim_)
: t(t_), imd(imd_), pai(pai_), U(U_), ctx(ctx_), qdim(qdim_),
- cv_old(-1) {}
+ cv_old(-1)
+ {}
};
struct ga_instruction_slice_local_dofs : public ga_instruction {
@@ -289,7 +290,7 @@ namespace getfem {
const mesh_im_level_set *mimls;
const fem_interpolation_context &ctx;
base_small_vector vec;
-
+
virtual int exec() {
GA_DEBUG_INFO("Instruction: unit normal vector to a level-set");
mimls->compute_normal_vector(ctx, vec);
@@ -1593,7 +1594,7 @@ namespace getfem {
return 0;
}
ga_instruction_add_to_coeff(base_tensor &t_, const base_tensor &tc1_,
- scalar_type &coeff_)
+ scalar_type &coeff_)
: t(t_), tc1(tc1_), coeff(coeff_) {}
};
@@ -1771,21 +1772,21 @@ namespace getfem {
size_type n0 = tc1.size() / (n1*n2*nn);
auto it = t.begin();
for (size_type i = 0; i < nn; ++i) {
- size_type s1 = i*n1*n2*n0;
- for (size_type j = 0; j < n1; ++j) {
- size_type s2 = s1 + j*n0;
- for (size_type k = 0; k < n2; ++k) {
- size_type s3 = s2 + k*n1*n0;
- for (size_type l = 0; l < n0; ++l, ++it)
- *it = tc1[s3+l];
- }
- }
+ size_type s1 = i*n1*n2*n0;
+ for (size_type j = 0; j < n1; ++j) {
+ size_type s2 = s1 + j*n0;
+ for (size_type k = 0; k < n2; ++k) {
+ size_type s3 = s2 + k*n1*n0;
+ for (size_type l = 0; l < n0; ++l, ++it)
+ *it = tc1[s3+l];
+ }
+ }
}
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_transpose(base_tensor &t_, const base_tensor &tc1_,
- size_type n1_, size_type n2_, size_type nn_)
+ size_type n1_, size_type n2_, size_type nn_)
: t(t_), tc1(tc1_), n1(n1_), n2(n2_), nn(nn_) {}
};
@@ -1800,19 +1801,19 @@ namespace getfem {
auto it = t.begin();
for (size_type i = 0; i < ii3; ++i)
- for (size_type j = 0; j < nn1; ++j)
- for (size_type k = 0; k < ii2; ++k)
- for (size_type l = 0; l < nn2; ++l) {
- size_type ind = j*ii1+k*ii1*nn1+l*ii1*nn1*ii2+i*ii1*nn1*ii2*nn2;
- for (size_type m = 0; m < ii1; ++m, ++it)
- *it = tc1[m+ind];
- }
+ for (size_type j = 0; j < nn1; ++j)
+ for (size_type k = 0; k < ii2; ++k)
+ for (size_type l = 0; l < nn2; ++l) {
+ size_type ind = j*ii1+k*ii1*nn1+l*ii1*nn1*ii2+i*ii1*nn1*ii2*nn2;
+ for (size_type m = 0; m < ii1; ++m, ++it)
+ *it = tc1[m+ind];
+ }
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_swap_indices(base_tensor &t_, const base_tensor &tc1_,
- size_type n1_, size_type n2_,
- size_type i2_, size_type i3_)
+ size_type n1_, size_type n2_,
+ size_type i2_, size_type i3_)
: t(t_), tc1(tc1_), nn1(n1_), nn2(n2_), ii2(i2_), ii3(i3_) {}
};
@@ -1827,19 +1828,19 @@ namespace getfem {
auto it = t.begin();
for (size_type i = 0; i < nn; ++i)
- for (size_type j = 0; j < ii2; ++j) {
- size_type ind = i*ii1+j*ii1*nn;
- for (size_type k = 0; k < ii1; ++k, ++it)
- *it = tc1[k+ind];
- }
+ for (size_type j = 0; j < ii2; ++j) {
+ size_type ind = i*ii1+j*ii1*nn;
+ for (size_type k = 0; k < ii1; ++k, ++it)
+ *it = tc1[k+ind];
+ }
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_index_move_last(base_tensor &t_, const base_tensor &tc1_,
- size_type n_, size_type i2_)
+ size_type n_, size_type i2_)
: t(t_), tc1(tc1_), nn(n_), ii2(i2_) {}
};
-
+
struct ga_instruction_transpose_no_test : public ga_instruction {
base_tensor &t;
const base_tensor &tc1;
@@ -1850,19 +1851,19 @@ namespace getfem {
auto it = t.begin();
for (size_type i = 0; i < nn; ++i) {
- size_type s1 = i*n1*n2;
- for (size_type j = 0; j < n1; ++j) {
- size_type s2 = s1 + j;
- for (size_type k = 0; k < n2; ++k, ++it)
- *it = tc1[s2 + k*n1];
- }
+ size_type s1 = i*n1*n2;
+ for (size_type j = 0; j < n1; ++j) {
+ size_type s2 = s1 + j;
+ for (size_type k = 0; k < n2; ++k, ++it)
+ *it = tc1[s2 + k*n1];
+ }
}
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_transpose_no_test(base_tensor &t_, const base_tensor &tc1_,
- size_type n1_, size_type n2_,
- size_type nn_)
+ size_type n1_, size_type n2_,
+ size_type nn_)
: t(t_), tc1(tc1_), n1(n1_), n2(n2_), nn(nn_) {}
};
@@ -2077,22 +2078,22 @@ namespace getfem {
GA_DEBUG_INFO("Instruction: single contraction on a single tensor");
size_type ii1 = tc1.size() / (nn*nn*ii2*ii3);
-
+
base_tensor::iterator it = t.begin();
for (size_type i = 0; i < ii3; ++i)
- for (size_type j = 0; j < ii2; ++j)
- for (size_type k = 0; k < ii1; ++k, ++it) {
- *it = scalar_type(0);
- size_type pre_ind = k+j*ii1*nn+i*ii1*nn*ii2*nn;
- for (size_type n = 0; n < nn; ++n)
- *it += tc1[pre_ind+n*ii1+n*ii1*nn*ii2];
- }
-
+ for (size_type j = 0; j < ii2; ++j)
+ for (size_type k = 0; k < ii1; ++k, ++it) {
+ *it = scalar_type(0);
+ size_type pre_ind = k+j*ii1*nn+i*ii1*nn*ii2*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += tc1[pre_ind+n*ii1+n*ii1*nn*ii2];
+ }
+
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_contract_1_1(base_tensor &t_, base_tensor &tc1_,
- size_type n_, size_type i2_, size_type i3_)
+ size_type n_, size_type i2_, size_type i3_)
: t(t_), tc1(tc1_), nn(n_), ii2(i2_), ii3(i3_) {}
};
@@ -2105,30 +2106,30 @@ namespace getfem {
size_type ift1 = tc1.size() / (nn*ii1*ii2);
size_type ift2 = tc2.size() / (nn*ii3*ii4);
-
+
base_tensor::iterator it = t.begin();
for (size_type i = 0; i < ii4; ++i)
- for (size_type j = 0; j < ii3; ++j)
- for (size_type k = 0; k < ii2; ++k)
- for (size_type l = 0; l < ii1; ++l)
- for (size_type p = 0; p < ift2; ++p)
- for (size_type q = 0; q < ift1; ++q, ++it) {
- *it = scalar_type(0);
- size_type ind1 = q+l*ift1+k*ift1*ii1*nn;
- size_type ind2 = p+j*ift2+i*ift2*ii3*nn;
- for (size_type n = 0; n < nn; ++n)
- *it += tc1[ind1+n*ift1*ii1] * tc2[ind2+n*ift2*ii3];
- }
-
+ for (size_type j = 0; j < ii3; ++j)
+ for (size_type k = 0; k < ii2; ++k)
+ for (size_type l = 0; l < ii1; ++l)
+ for (size_type p = 0; p < ift2; ++p)
+ for (size_type q = 0; q < ift1; ++q, ++it) {
+ *it = scalar_type(0);
+ size_type ind1 = q+l*ift1+k*ift1*ii1*nn;
+ size_type ind2 = p+j*ift2+i*ift2*ii3*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += tc1[ind1+n*ift1*ii1] * tc2[ind2+n*ift2*ii3];
+ }
+
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_contract_2_1(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_,
- size_type n_, size_type i1_, size_type i2_,
- size_type i3_, size_type i4_)
+ base_tensor &tc2_,
+ size_type n_, size_type i1_, size_type i2_,
+ size_type i3_, size_type i4_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_),
- ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_) {}
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_) {}
};
// Performs Amijk Bnljp -> Cnmiklp. To be optimized
@@ -2140,30 +2141,30 @@ namespace getfem {
size_type ift1 = tc1.size() / (nn*ii1*ii2);
size_type ift2 = tc2.size() / (nn*ii3*ii4);
-
+
base_tensor::iterator it = t.begin();
for (size_type i = 0; i < ii4; ++i)
- for (size_type j = 0; j < ii3; ++j)
- for (size_type k = 0; k < ii2; ++k)
- for (size_type l = 0; l < ii1; ++l)
- for (size_type q = 0; q < ift1; ++q)
- for (size_type p = 0; p < ift2; ++p, ++it) {
- *it = scalar_type(0);
- size_type ind1 = q+l*ift1+k*ift1*ii1*nn;
- size_type ind2 = p+j*ift2+i*ift2*ii3*nn;
- for (size_type n = 0; n < nn; ++n)
- *it += tc1[ind1+n*ift1*ii1] * tc2[ind2+n*ift2*ii3];
- }
-
+ for (size_type j = 0; j < ii3; ++j)
+ for (size_type k = 0; k < ii2; ++k)
+ for (size_type l = 0; l < ii1; ++l)
+ for (size_type q = 0; q < ift1; ++q)
+ for (size_type p = 0; p < ift2; ++p, ++it) {
+ *it = scalar_type(0);
+ size_type ind1 = q+l*ift1+k*ift1*ii1*nn;
+ size_type ind2 = p+j*ift2+i*ift2*ii3*nn;
+ for (size_type n = 0; n < nn; ++n)
+ *it += tc1[ind1+n*ift1*ii1] * tc2[ind2+n*ift2*ii3];
+ }
+
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_contract_2_1_rev(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_,
- size_type n_, size_type i1_, size_type i2_,
- size_type i3_, size_type i4_)
+ base_tensor &tc2_,
+ size_type n_, size_type i1_, size_type i2_,
+ size_type i3_, size_type i4_)
: t(t_), tc1(tc1_), tc2(tc2_), nn(n_),
- ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_) {}
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_) {}
};
// Performs Amijklp Bnqjrls -> Cmnikpqrs. To be optimized
@@ -2179,39 +2180,39 @@ namespace getfem {
size_type sn1 = ift2*ii4, sn2 = ift2*ii4*nn1*ii5;
if (inv_tc2) std::swap(sn1, sn2);
-
+
base_tensor::iterator it = t.begin();
for (size_type i = 0; i < ii6; ++i)
- for (size_type j = 0; j < ii5; ++j)
- for (size_type k = 0; k < ii4; ++k)
- for (size_type l = 0; l < ii3; ++l)
- for (size_type p = 0; p < ii2; ++p)
- for (size_type q = 0; q < ii1; ++q)
- for (size_type r = 0; r < ift2; ++r)
- for (size_type s = 0; s < ift1; ++s, ++it) {
- *it = scalar_type(0);
- size_type ind1
- = s+q*ift1+p*ift1*ii1*nn1+l*ift1*ii1*nn1*ii2*nn2;
- size_type ind2
- = r+k*ift2+j*ift2*ii4*nn1+i*ift2*ii4*nn1*ii5*nn2;
- for (size_type n1 = 0; n1 < nn1; ++n1)
- for (size_type n2 = 0; n2 < nn2; ++n2)
- *it += tc1[ind1+n1*ift1*ii1+n2*ift1*ii1*nn1*ii2]
- * tc2[ind2+n1*sn1+n2*sn2];
- }
-
+ for (size_type j = 0; j < ii5; ++j)
+ for (size_type k = 0; k < ii4; ++k)
+ for (size_type l = 0; l < ii3; ++l)
+ for (size_type p = 0; p < ii2; ++p)
+ for (size_type q = 0; q < ii1; ++q)
+ for (size_type r = 0; r < ift2; ++r)
+ for (size_type s = 0; s < ift1; ++s, ++it) {
+ *it = scalar_type(0);
+ size_type ind1
+ = s+q*ift1+p*ift1*ii1*nn1+l*ift1*ii1*nn1*ii2*nn2;
+ size_type ind2
+ = r+k*ift2+j*ift2*ii4*nn1+i*ift2*ii4*nn1*ii5*nn2;
+ for (size_type n1 = 0; n1 < nn1; ++n1)
+ for (size_type n2 = 0; n2 < nn2; ++n2)
+ *it += tc1[ind1+n1*ift1*ii1+n2*ift1*ii1*nn1*ii2]
+ * tc2[ind2+n1*sn1+n2*sn2];
+ }
+
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_contract_2_2(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_,
- size_type n1_, size_type n2_,
- size_type i1_, size_type i2_, size_type i3_,
- size_type i4_, size_type i5_, size_type i6_,
- bool intc2)
+ base_tensor &tc2_,
+ size_type n1_, size_type n2_,
+ size_type i1_, size_type i2_, size_type i3_,
+ size_type i4_, size_type i5_, size_type i6_,
+ bool intc2)
: t(t_), tc1(tc1_), tc2(tc2_), nn1(n1_), nn2(n2_),
- ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_), ii5(i5_), ii6(i6_),
- inv_tc2(intc2) {}
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_), ii5(i5_), ii6(i6_),
+ inv_tc2(intc2) {}
};
// Performs Amijklp Bnqjrls -> Cnmikpqrs. To be optimized
@@ -2227,39 +2228,39 @@ namespace getfem {
size_type sn1 = ift2*ii4, sn2 = ift2*ii4*nn1*ii5;
if (inv_tc2) std::swap(sn1, sn2);
-
+
base_tensor::iterator it = t.begin();
for (size_type i = 0; i < ii6; ++i)
- for (size_type j = 0; j < ii5; ++j)
- for (size_type k = 0; k < ii4; ++k)
- for (size_type l = 0; l < ii3; ++l)
- for (size_type p = 0; p < ii2; ++p)
- for (size_type q = 0; q < ii1; ++q)
- for (size_type s = 0; s < ift1; ++s)
- for (size_type r = 0; r < ift2; ++r, ++it) {
- *it = scalar_type(0);
- size_type ind1
- = s+q*ift1+p*ift1*ii1*nn1+l*ift1*ii1*nn1*ii2*nn2;
- size_type ind2
- = r+k*ift2+j*ift2*ii4*nn1+i*ift2*ii4*nn1*ii5*nn2;
- for (size_type n1 = 0; n1 < nn1; ++n1)
- for (size_type n2 = 0; n2 < nn2; ++n2)
- *it += tc1[ind1+n1*ift1*ii1+n2*ift1*ii1*nn1*ii2]
- * tc2[ind2+n1*sn1+n2*sn2];
- }
-
+ for (size_type j = 0; j < ii5; ++j)
+ for (size_type k = 0; k < ii4; ++k)
+ for (size_type l = 0; l < ii3; ++l)
+ for (size_type p = 0; p < ii2; ++p)
+ for (size_type q = 0; q < ii1; ++q)
+ for (size_type s = 0; s < ift1; ++s)
+ for (size_type r = 0; r < ift2; ++r, ++it) {
+ *it = scalar_type(0);
+ size_type ind1
+ = s+q*ift1+p*ift1*ii1*nn1+l*ift1*ii1*nn1*ii2*nn2;
+ size_type ind2
+ = r+k*ift2+j*ift2*ii4*nn1+i*ift2*ii4*nn1*ii5*nn2;
+ for (size_type n1 = 0; n1 < nn1; ++n1)
+ for (size_type n2 = 0; n2 < nn2; ++n2)
+ *it += tc1[ind1+n1*ift1*ii1+n2*ift1*ii1*nn1*ii2]
+ * tc2[ind2+n1*sn1+n2*sn2];
+ }
+
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_contract_2_2_rev(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_,
- size_type n1_, size_type n2_,
- size_type i1_, size_type i2_, size_type i3_,
- size_type i4_, size_type i5_, size_type i6_,
- bool intc2)
+ base_tensor &tc2_,
+ size_type n1_, size_type n2_,
+ size_type i1_, size_type i2_, size_type
i3_,
+ size_type i4_, size_type i5_, size_type
i6_,
+ bool intc2)
: t(t_), tc1(tc1_), tc2(tc2_), nn1(n1_), nn2(n2_),
- ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_), ii5(i5_), ii6(i6_),
- inv_tc2(intc2) {}
+ ii1(i1_), ii2(i2_), ii3(i3_), ii4(i4_), ii5(i5_), ii6(i6_),
+ inv_tc2(intc2) {}
};
@@ -2269,7 +2270,7 @@ namespace getfem {
size_type n;
virtual int exec() {
GA_DEBUG_INFO("Instruction: order one contraction "
- "(dot product or matrix multiplication)");
+ "(dot product or matrix multiplication)");
size_type s1 = tc1.size() / n;
size_type s2 = tc2.size() / n;
@@ -2296,25 +2297,25 @@ namespace getfem {
// t of size q*l*m*p
virtual int exec() {
GA_DEBUG_INFO("Instruction: specific order one contraction "
- "(dot product or matrix multiplication)");
+ "(dot product or matrix multiplication)");
size_type q = tc1.size() / (m * n);
size_type l = tc2.size() / (p * n);
-
+
base_tensor::iterator it = t.begin();
for (size_type r = 0; r < p; ++r)
- for (size_type k = 0; k < m; ++k)
- for (size_type j = 0; j < l; ++j)
- for (size_type i = 0; i < q; ++i, ++it) {
- *it = scalar_type(0);
- for (size_type s = 0; s < n; ++s)
- *it += tc1[i+k*q+s*q*m] * tc2[j+s*l+r*l*n];
- }
+ for (size_type k = 0; k < m; ++k)
+ for (size_type j = 0; j < l; ++j)
+ for (size_type i = 0; i < q; ++i, ++it) {
+ *it = scalar_type(0);
+ for (size_type s = 0; s < n; ++s)
+ *it += tc1[i+k*q+s*q*m] * tc2[j+s*l+r*l*n];
+ }
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_matrix_mult_spec(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_, size_type n_,
- size_type m_, size_type p_)
+ base_tensor &tc2_, size_type n_,
+ size_type m_, size_type p_)
: t(t_), tc1(tc1_), tc2(tc2_), n(n_), m(m_), p(p_) {}
};
@@ -2325,25 +2326,25 @@ namespace getfem {
// t of size l*q*m*p
virtual int exec() {
GA_DEBUG_INFO("Instruction: specific order one contraction "
- "(dot product or matrix multiplication)");
+ "(dot product or matrix multiplication)");
size_type q = tc1.size() / (m * n);
size_type l = tc2.size() / (p * n);
-
+
base_tensor::iterator it = t.begin();
for (size_type r = 0; r < p; ++r)
- for (size_type k = 0; k < m; ++k)
- for (size_type i = 0; i < q; ++i)
- for (size_type j = 0; j < l; ++j, ++it) {
- *it = scalar_type(0);
- for (size_type s = 0; s < n; ++s)
- *it += tc1[i+k*q+s*q*m] * tc2[j+s*l+r*l*n];
- }
+ for (size_type k = 0; k < m; ++k)
+ for (size_type i = 0; i < q; ++i)
+ for (size_type j = 0; j < l; ++j, ++it) {
+ *it = scalar_type(0);
+ for (size_type s = 0; s < n; ++s)
+ *it += tc1[i+k*q+s*q*m] * tc2[j+s*l+r*l*n];
+ }
GA_DEBUG_ASSERT(it == t.end(), "Wrong sizes");
return 0;
}
ga_instruction_matrix_mult_spec2(base_tensor &t_, base_tensor &tc1_,
- base_tensor &tc2_, size_type n_,
- size_type m_, size_type p_)
+ base_tensor &tc2_, size_type n_,
+ size_type m_, size_type p_)
: t(t_), tc1(tc1_), tc2(tc2_), n(n_), m(m_), p(p_) {}
};
@@ -4675,14 +4676,14 @@ namespace getfem {
add_interval_to_gis(workspace, v, gis);
} else {
if (gis.var_intervals.find(varname) == gis.var_intervals.end()) {
- const mesh_fem *mf = workspace.associated_mf(varname);
- size_type nd = mf ? mf->nb_basic_dof() :
- gmm::vect_size(workspace.value(varname));
- gis.var_intervals[varname]=gmm::sub_interval(gis.nb_dof, nd);
- gis.nb_dof += nd;
+ const mesh_fem *mf = workspace.associated_mf(varname);
+ size_type nd = mf ? mf->nb_basic_dof() :
+ gmm::vect_size(workspace.value(varname));
+ gis.var_intervals[varname]=gmm::sub_interval(gis.nb_dof, nd);
+ gis.nb_dof += nd;
}
gis.max_dof = std::max(gis.max_dof,
- workspace.interval_of_variable(varname).last());
+ workspace.interval_of_variable(varname).last());
}
}
@@ -4955,22 +4956,22 @@ namespace getfem {
case GA_NODE_NORMAL:
{
- GMM_ASSERT1(!function_case,
- "No use of Normal is allowed in functions");
- if (pnode->tensor().size() != m.dim())
- pnode->init_vector_tensor(m.dim());
- const mesh_im_level_set *mimls
- = dynamic_cast<const mesh_im_level_set *>(rmi.im);
- if (mimls && mimls->location()==mesh_im_level_set::INTEGRATE_BOUNDARY) {
- // Appel avec ctx (pt de Gauss)
- pgai = std::make_shared<ga_instruction_level_set_normal_vector>
- (pnode->tensor(), mimls, gis.ctx);
- rmi.instructions.push_back(std::move(pgai));
- } else {
- pgai = std::make_shared<ga_instruction_copy_Normal>
- (pnode->tensor(), gis.Normal);
- rmi.instructions.push_back(std::move(pgai));
- }
+ GMM_ASSERT1(!function_case,
+ "No use of Normal is allowed in functions");
+ if (pnode->tensor().size() != m.dim())
+ pnode->init_vector_tensor(m.dim());
+ const mesh_im_level_set *mimls
+ = dynamic_cast<const mesh_im_level_set *>(rmi.im);
+ if (mimls && mimls->location()==mesh_im_level_set::INTEGRATE_BOUNDARY)
{
+ // Appel avec ctx (pt de Gauss)
+ pgai = std::make_shared<ga_instruction_level_set_normal_vector>
+ (pnode->tensor(), mimls, gis.ctx);
+ rmi.instructions.push_back(std::move(pgai));
+ } else {
+ pgai = std::make_shared<ga_instruction_copy_Normal>
+ (pnode->tensor(), gis.Normal);
+ rmi.instructions.push_back(std::move(pgai));
+ }
}
break;
@@ -5681,7 +5682,7 @@ namespace getfem {
pgai = pga_instruction();
if ((pnode->op_type == GA_DOT && dim1 <= 1) ||
- pnode->op_type == GA_COLON ||
+ pnode->op_type == GA_COLON ||
(pnode->op_type == GA_MULT && dim0 == 4) ||
(pnode->op_type == GA_MULT && dim1 <= 1) ||
child0->tensor().size() == 1 || tps1 == 1) {
@@ -5775,7 +5776,7 @@ namespace getfem {
}
}
} else { // GA_MULT or GA_DOT for dim1 > 1
- // and child1->tensor_proper_size() > 1
+ // and child1->tensor_proper_size() > 1
if (pnode->test_function_type < 3) {
if (tps0 == 1) {
if (is_uniform) // Unrolled instruction
@@ -5785,13 +5786,13 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_simple_tmult>
(pnode->tensor(), child0->tensor(), child1->tensor());
} else {
- if (child1->test_function_type == 0)
- pgai = std::make_shared<ga_instruction_matrix_mult>
- (pnode->tensor(), child0->tensor(), child1->tensor(), s2);
- else
- pgai = std::make_shared<ga_instruction_matrix_mult_spec>
- (pnode->tensor(), child0->tensor(), child1->tensor(),
- s2, tps0/s2, tps1/s2);
+ if (child1->test_function_type == 0)
+ pgai = std::make_shared<ga_instruction_matrix_mult>
+ (pnode->tensor(), child0->tensor(), child1->tensor(), s2);
+ else
+ pgai = std::make_shared<ga_instruction_matrix_mult_spec>
+ (pnode->tensor(), child0->tensor(), child1->tensor(),
+ s2, tps0/s2, tps1/s2);
}
} else {
if (child0->tensor_proper_size() == 1) {
@@ -5812,13 +5813,13 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_matrix_mult>
(pnode->tensor(), child0->tensor(), child1->tensor(), s2);
else if (child1->test_function_type == 2)
- pgai = std::make_shared<ga_instruction_matrix_mult_spec>
+ pgai = std::make_shared<ga_instruction_matrix_mult_spec>
(pnode->tensor(), child0->tensor(), child1->tensor(),
- s2, tps0/s2, tps1/s2);
- else
- pgai = std::make_shared<ga_instruction_matrix_mult_spec2>
+ s2, tps0/s2, tps1/s2);
+ else
+ pgai = std::make_shared<ga_instruction_matrix_mult_spec2>
(pnode->tensor(), child0->tensor(), child1->tensor(),
- s2, tps0/s2, tps1/s2);
+ s2, tps0/s2, tps1/s2);
}
}
}
@@ -5847,18 +5848,18 @@ namespace getfem {
case GA_QUOTE:
if (pnode->tensor_proper_size() > 1) {
- size_type n1 = child0->tensor_proper_size(0);
- size_type n2 = (child0->tensor_order() > 1) ?
- child0->tensor_proper_size(1) : 1;
- size_type nn = 1;
- for (size_type i = 2; i < child0->tensor_order(); ++i)
- nn *= child0->tensor_proper_size(i);
- if (child0->nb_test_functions() == 0)
- pgai = std::make_shared<ga_instruction_transpose_no_test>
- (pnode->tensor(), child0->tensor(), n1, n2, nn);
- else
- pgai = std::make_shared<ga_instruction_transpose>
- (pnode->tensor(), child0->tensor(), n1, n2, nn);
+ size_type n1 = child0->tensor_proper_size(0);
+ size_type n2 = (child0->tensor_order() > 1) ?
+ child0->tensor_proper_size(1) : 1;
+ size_type nn = 1;
+ for (size_type i = 2; i < child0->tensor_order(); ++i)
+ nn *= child0->tensor_proper_size(i);
+ if (child0->nb_test_functions() == 0)
+ pgai = std::make_shared<ga_instruction_transpose_no_test>
+ (pnode->tensor(), child0->tensor(), n1, n2, nn);
+ else
+ pgai = std::make_shared<ga_instruction_transpose>
+ (pnode->tensor(), child0->tensor(), n1, n2, nn);
rmi.instructions.push_back(std::move(pgai));
} else {
pnode->t.set_to_copy(child0->t);
@@ -5994,14 +5995,14 @@ namespace getfem {
{
if (pnode->test_function_type) {
std::vector<const base_tensor *> components(pnode->children.size());
- for (size_type i = 0; i < pnode->children.size(); ++i)
- components[i] = &(pnode->children[i]->tensor());
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ components[i] = &(pnode->children[i]->tensor());
pgai = std::make_shared<ga_instruction_c_matrix_with_tests>
(pnode->tensor(), components);
} else {
std::vector<scalar_type *> components(pnode->children.size());
- for (size_type i = 0; i < pnode->children.size(); ++i)
- components[i] = &(pnode->children[i]->tensor()[0]);
+ for (size_type i = 0; i < pnode->children.size(); ++i)
+ components[i] = &(pnode->children[i]->tensor()[0]);
pgai = std::make_shared<ga_instruction_simple_c_matrix>
(pnode->tensor(), components);
}
@@ -6015,112 +6016,112 @@ namespace getfem {
child1->tensor());
rmi.instructions.push_back(std::move(pgai));
} else if (child0->node_type == GA_NODE_IND_MOVE_LAST) {
- size_type ind;
- ind = size_type(round(pnode->children[2]->tensor()[0])-1);
- size_type ii2 = 1;
- for (size_type i = 0; i < child1->tensor_order(); ++i)
- if (i>ind) ii2 *= child1->tensor_proper_size(i);
- size_type nn = child1->tensor_proper_size(ind);
+ size_type ind;
+ ind = size_type(round(pnode->children[2]->tensor()[0])-1);
+ size_type ii2 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i)
+ if (i>ind) ii2 *= child1->tensor_proper_size(i);
+ size_type nn = child1->tensor_proper_size(ind);
pgai = std::make_shared<ga_instruction_index_move_last>
- (pnode->tensor(), child1->tensor(), nn, ii2);
+ (pnode->tensor(), child1->tensor(), nn, ii2);
rmi.instructions.push_back(std::move(pgai));
} else if (child0->node_type == GA_NODE_SWAP_IND) {
- size_type ind[4];
- for (size_type i = 2; i < 4; ++i)
- ind[i] = size_type(round(pnode->children[i]->tensor()[0])-1);
- if (ind[2] > ind[3]) std::swap(ind[2], ind[3]);
- size_type ii2 = 1, ii3 = 1;
- for (size_type i = 0; i < child1->tensor_order(); ++i) {
- if (i>ind[2] && i<ind[3]) ii2 *= child1->tensor_proper_size(i);
- if (i>ind[3]) ii3 *= child1->tensor_proper_size(i);
- }
- size_type nn1 = child1->tensor_proper_size(ind[2]);
- size_type nn2 = child1->tensor_proper_size(ind[3]);
-
+ size_type ind[4];
+ for (size_type i = 2; i < 4; ++i)
+ ind[i] = size_type(round(pnode->children[i]->tensor()[0])-1);
+ if (ind[2] > ind[3]) std::swap(ind[2], ind[3]);
+ size_type ii2 = 1, ii3 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i>ind[2] && i<ind[3]) ii2 *= child1->tensor_proper_size(i);
+ if (i>ind[3]) ii3 *= child1->tensor_proper_size(i);
+ }
+ size_type nn1 = child1->tensor_proper_size(ind[2]);
+ size_type nn2 = child1->tensor_proper_size(ind[3]);
+
pgai = std::make_shared<ga_instruction_swap_indices>
- (pnode->tensor(), child1->tensor(), nn1, nn2, ii2, ii3);
+ (pnode->tensor(), child1->tensor(), nn1, nn2, ii2, ii3);
rmi.instructions.push_back(std::move(pgai));
} else if (child0->node_type == GA_NODE_CONTRACT) {
- std::vector<size_type> ind(2), indsize(2);
- pga_tree_node child2(0);
- if (pnode->children.size() == 4)
- { ind[0] = 2; ind[1] = 3; }
- else if (pnode->children.size() == 5)
- { ind[0] = 2; ind[1] = 4; child2 = pnode->children[3]; }
- else if (pnode->children.size() == 7) {
- ind.resize(4); indsize.resize(4);
- ind[0] = 2; ind[1] = 3; ind[2] = 5; ind[3] = 6;
- child2 = pnode->children[4];
- }
- size_type kk = 0, ll = 1;
- for (size_type i = 1; i < pnode->children.size(); ++i) {
- if (i == ind[kk]) {
- ind[kk] = size_type(round(pnode->children[i]->tensor()[0])-1);
- indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
- ++kk;
- } else ll = i;
- }
-
- if (pnode->children.size() == 4) {
- size_type i1 = ind[0], i2 = ind[1];
- if (i1 > i2) std::swap(i1, i2);
- size_type ii2 = 1, ii3 = 1;
- for (size_type i = 0; i < child1->tensor_order(); ++i) {
- if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
- if (i > i2) ii3 *= child1->tensor_proper_size(i);
- }
- pgai = std::make_shared<ga_instruction_contract_1_1>
- (pnode->tensor(), child1->tensor(), indsize[0], ii2, ii3);
- }
- else if (pnode->children.size() == 5) {
- // Particular cases should be detected (ii2=ii3=1 in particular).
- size_type i1 = ind[0], i2 = ind[1];
- size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1;
- for (size_type i = 0; i < child1->tensor_order(); ++i) {
- if (i < i1) ii1 *= child1->tensor_proper_size(i);
- if (i > i1) ii2 *= child1->tensor_proper_size(i);
- }
- for (size_type i = 0; i < child2->tensor_order(); ++i) {
- if (i < i2) ii3 *= child2->tensor_proper_size(i);
- if (i > i2) ii4 *= child2->tensor_proper_size(i);
- }
- if (child1->test_function_type==1 && child2->test_function_type==2)
- pgai = std::make_shared<ga_instruction_contract_2_1_rev>
- (pnode->tensor(), child1->tensor(), child2->tensor(),
- indsize[0], ii1, ii2, ii3, ii4);
- else
- pgai = std::make_shared<ga_instruction_contract_2_1>
- (pnode->tensor(), child1->tensor(), child2->tensor(),
- indsize[0], ii1, ii2, ii3, ii4);
- }
- else if (pnode->children.size() == 7) {
- // Particular cases should be detected (ii2=ii3=1 in particular).
- size_type i1 = ind[0], i2 = ind[1], i3 = ind[2], i4 = ind[3];
- size_type nn1 = indsize[0], nn2 = indsize[1];
- size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1, ii5 = 1, ii6 = 1;
- if (i1 > i2)
- { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
- for (size_type i = 0; i < child1->tensor_order(); ++i) {
- if (i < i1) ii1 *= child1->tensor_proper_size(i);
- if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
- if (i > i2) ii3 *= child1->tensor_proper_size(i);
- }
- for (size_type i = 0; i < child2->tensor_order(); ++i) {
- if (i < i3 && i < i4) ii4 *= child2->tensor_proper_size(i);
- if ((i > i3 && i < i4) || (i > i4 && i < i3))
- ii5 *= child2->tensor_proper_size(i);
- if (i > i3 && i > i4) ii6 *= child2->tensor_proper_size(i);
- }
- if (child1->test_function_type==1 && child2->test_function_type==2)
- pgai = std::make_shared<ga_instruction_contract_2_2_rev>
- (pnode->tensor(), child1->tensor(), child2->tensor(),
- nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6, i4 < i3);
- else
- pgai = std::make_shared<ga_instruction_contract_2_2>
- (pnode->tensor(), child1->tensor(), child2->tensor(),
- nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6, i4 < i3);
- }
- rmi.instructions.push_back(std::move(pgai));
+ std::vector<size_type> ind(2), indsize(2);
+ pga_tree_node child2(0);
+ if (pnode->children.size() == 4)
+ { ind[0] = 2; ind[1] = 3; }
+ else if (pnode->children.size() == 5)
+ { ind[0] = 2; ind[1] = 4; child2 = pnode->children[3]; }
+ else if (pnode->children.size() == 7) {
+ ind.resize(4); indsize.resize(4);
+ ind[0] = 2; ind[1] = 3; ind[2] = 5; ind[3] = 6;
+ child2 = pnode->children[4];
+ }
+ size_type kk = 0, ll = 1;
+ for (size_type i = 1; i < pnode->children.size(); ++i) {
+ if (i == ind[kk]) {
+ ind[kk] = size_type(round(pnode->children[i]->tensor()[0])-1);
+ indsize[kk] = pnode->children[ll]->tensor_proper_size(ind[kk]);
+ ++kk;
+ } else ll = i;
+ }
+
+ if (pnode->children.size() == 4) {
+ size_type i1 = ind[0], i2 = ind[1];
+ if (i1 > i2) std::swap(i1, i2);
+ size_type ii2 = 1, ii3 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
+ if (i > i2) ii3 *= child1->tensor_proper_size(i);
+ }
+ pgai = std::make_shared<ga_instruction_contract_1_1>
+ (pnode->tensor(), child1->tensor(), indsize[0], ii2, ii3);
+ }
+ else if (pnode->children.size() == 5) {
+ // Particular cases should be detected (ii2=ii3=1 in particular).
+ size_type i1 = ind[0], i2 = ind[1];
+ size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1;
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1) ii2 *= child1->tensor_proper_size(i);
+ }
+ for (size_type i = 0; i < child2->tensor_order(); ++i) {
+ if (i < i2) ii3 *= child2->tensor_proper_size(i);
+ if (i > i2) ii4 *= child2->tensor_proper_size(i);
+ }
+ if (child1->test_function_type==1 && child2->test_function_type==2)
+ pgai = std::make_shared<ga_instruction_contract_2_1_rev>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ indsize[0], ii1, ii2, ii3, ii4);
+ else
+ pgai = std::make_shared<ga_instruction_contract_2_1>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ indsize[0], ii1, ii2, ii3, ii4);
+ }
+ else if (pnode->children.size() == 7) {
+ // Particular cases should be detected (ii2=ii3=1 in particular).
+ size_type i1 = ind[0], i2 = ind[1], i3 = ind[2], i4 = ind[3];
+ size_type nn1 = indsize[0], nn2 = indsize[1];
+ size_type ii1 = 1, ii2 = 1, ii3 = 1, ii4 = 1, ii5 = 1, ii6 = 1;
+ if (i1 > i2)
+ { std::swap(i1, i2); std::swap(i3, i4); std::swap(nn1, nn2); }
+ for (size_type i = 0; i < child1->tensor_order(); ++i) {
+ if (i < i1) ii1 *= child1->tensor_proper_size(i);
+ if (i > i1 && i < i2) ii2 *= child1->tensor_proper_size(i);
+ if (i > i2) ii3 *= child1->tensor_proper_size(i);
+ }
+ for (size_type i = 0; i < child2->tensor_order(); ++i) {
+ if (i < i3 && i < i4) ii4 *= child2->tensor_proper_size(i);
+ if ((i > i3 && i < i4) || (i > i4 && i < i3))
+ ii5 *= child2->tensor_proper_size(i);
+ if (i > i3 && i > i4) ii6 *= child2->tensor_proper_size(i);
+ }
+ if (child1->test_function_type==1 && child2->test_function_type==2)
+ pgai = std::make_shared<ga_instruction_contract_2_2_rev>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6, i4 < i3);
+ else
+ pgai = std::make_shared<ga_instruction_contract_2_2>
+ (pnode->tensor(), child1->tensor(), child2->tensor(),
+ nn1, nn2, ii1, ii2, ii3, ii4, ii5, ii6, i4 < i3);
+ }
+ rmi.instructions.push_back(std::move(pgai));
} else if (child0->node_type == GA_NODE_PREDEF_FUNC) {
std::string name = child0->name;
@@ -6252,7 +6253,7 @@ namespace getfem {
}
void ga_compile_function(ga_workspace &workspace,
- ga_instruction_set &gis, bool scalar) {
+ ga_instruction_set &gis, bool scalar) {
for (size_type i = 0; i < workspace.nb_trees(); ++i) {
const ga_workspace::tree_description &td = workspace.tree_info(i);
@@ -6269,9 +6270,9 @@ namespace getfem {
gis.coeff = scalar_type(1);
pga_instruction pgai;
- workspace.assembled_tensor() = root->tensor();
- pgai = std::make_shared<ga_instruction_add_to_coeff>
- (workspace.assembled_tensor(), root->tensor(), gis.coeff);
+ workspace.assembled_tensor() = root->tensor();
+ pgai = std::make_shared<ga_instruction_add_to_coeff>
+ (workspace.assembled_tensor(), root->tensor(), gis.coeff);
gis.whole_instructions[rm].instructions.push_back(std::move(pgai));
}
}
@@ -6364,7 +6365,7 @@ namespace getfem {
}
void ga_compile_interpolation(ga_workspace &workspace,
- ga_instruction_set &gis) {
+ ga_instruction_set &gis) {
gis.transformations.clear();
gis.whole_instructions.clear();
for (size_type i = 0; i < workspace.nb_trees(); ++i) {
@@ -6383,7 +6384,8 @@ namespace getfem {
ga_instruction_set::region_mim rm(td.mim, td.rg);
ga_instruction_set::region_mim_instructions &rmi
= gis.whole_instructions[rm];
- rmi.m = td.m; rmi.im = td.mim;
+ rmi.m = td.m;
+ rmi.im = td.mim;
// rmi.interpolate_infos.clear();
ga_compile_interpolate_trans(root, workspace, gis, rmi, *(td.m));
ga_compile_node(root, workspace, gis,rmi, *(td.m), false,
@@ -6400,7 +6402,7 @@ namespace getfem {
}
void ga_compile(ga_workspace &workspace,
- ga_instruction_set &gis, size_type order) {
+ ga_instruction_set &gis, size_type order) {
gis.transformations.clear();
gis.whole_instructions.clear();
for (size_type version : std::array<size_type, 3>{1, 0, 2}) {
@@ -6433,7 +6435,8 @@ namespace getfem {
ga_instruction_set::region_mim rm(td.mim, td.rg);
ga_instruction_set::region_mim_instructions &rmi
= gis.whole_instructions[rm];
- rmi.m = td.m; rmi.im = td.mim;
+ rmi.m = td.m;
+ rmi.im = td.mim;
// rmi.interpolate_infos.clear();
ga_compile_interpolate_trans(root, workspace, gis, rmi, *(td.m));
ga_compile_node(root, workspace, gis, rmi, *(td.m), false,
@@ -6446,25 +6449,25 @@ namespace getfem {
auto *imd
= workspace.associated_im_data(td.varname_interpolation);
auto &V = const_cast<model_real_plain_vector &>
- (workspace.value(td.varname_interpolation));
+ (workspace.value(td.varname_interpolation));
GMM_ASSERT1(imd, "Internal error");
auto pgai = std::make_shared<ga_instruction_assignment>
- (root->tensor(), V, gis.ctx, imd);
+ (root->tensor(), V, gis.ctx, imd);
rmi.instructions.push_back(std::move(pgai));
- }
+ }
} else { // Addition of an assembly instruction
pga_instruction pgai;
switch(order) {
case 0:
- workspace.assembled_tensor() = root->tensor();
- pgai = std::make_shared<ga_instruction_add_to_coeff>
- (workspace.assembled_tensor(), root->tensor(), gis.coeff);
+ workspace.assembled_tensor() = root->tensor();
+ pgai = std::make_shared<ga_instruction_add_to_coeff>
+ (workspace.assembled_tensor(), root->tensor(), gis.coeff);
break;
case 1:
{
- GMM_ASSERT1(root->tensor_proper_size() == 1,
- "Invalid vector or tensor quantity. An order 1 "
- "weak form has to be a scalar quantity");
+ GMM_ASSERT1(root->tensor_proper_size() == 1,
+ "Invalid vector or tensor quantity. An order 1 "
+ "weak form has to be a scalar quantity");
const mesh_fem *mf=workspace.associated_mf(root->name_test1);
const mesh_fem **mfg = 0;
add_interval_to_gis(workspace, root->name_test1, gis);
@@ -6504,9 +6507,9 @@ namespace getfem {
break;
case 2:
{
- GMM_ASSERT1(root->tensor_proper_size() == 1,
- "Invalid vector or tensor quantity. An order 2 "
- "weak form has to be a scalar quantity");
+ GMM_ASSERT1(root->tensor_proper_size() == 1,
+ "Invalid vector or tensor quantity. An order 2 "
+ "weak form has to be a scalar quantity");
const mesh_fem
*mf1=workspace.associated_mf(root->name_test1);
const mesh_fem
*mf2=workspace.associated_mf(root->name_test2);
const mesh_fem **mfg1 = 0, **mfg2 = 0;
@@ -6627,8 +6630,8 @@ namespace getfem {
}
void ga_interpolation_exec(ga_instruction_set &gis,
- ga_workspace &workspace,
- ga_interpolation_context &gic) {
+ ga_workspace &workspace,
+ ga_interpolation_context &gic) {
base_matrix G;
base_small_vector un, up;
@@ -6824,7 +6827,7 @@ namespace getfem {
first_ind = pai->ind_first_point_on_face(v.f());
}
for (gis.ipt = 0; gis.ipt < gis.nbpt; ++(gis.ipt)) {
- // cout << "Gauss pt " << gis.ipt << endl;
+ // cout << "Gauss pt " << gis.ipt << endl;
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())) {
@@ -6841,11 +6844,11 @@ namespace getfem {
gmm::clean(gis.Normal, 1e-13);
} else gis.Normal.resize(0);
}
- auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
+ auto ipt_coeff = pai->coeff(first_ind+gis.ipt);
gis.coeff = J * ipt_coeff;
- bool enable_ipt = (gmm::abs(ipt_coeff) > 0.0 ||
- workspace.include_empty_int_points());
- if (!enable_ipt) gis.coeff = scalar_type(0);
+ 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;
diff --git a/src/getfem_generic_assembly_interpolation.cc
b/src/getfem_generic_assembly_interpolation.cc
index fcff8f7..2a48225 100644
--- a/src/getfem_generic_assembly_interpolation.cc
+++ b/src/getfem_generic_assembly_interpolation.cc
@@ -479,7 +479,7 @@ namespace getfem {
// Interpolate transformation with an expression
//=========================================================================
- class interpolate_transformation_expression
+ class interpolate_transformation_expression
: public virtual_interpolate_transformation, public context_dependencies {
struct workspace_gis_pair : public std::pair<ga_workspace,
ga_instruction_set> {
@@ -620,7 +620,7 @@ namespace getfem {
local_gis = ga_instruction_set();
}
- std::string expression(void) const { return expr; }
+ std::string expression() const { return expr; }
int transform(const ga_workspace &/*workspace*/, const mesh &m,
fem_interpolation_context &ctx_x,
@@ -637,7 +637,7 @@ namespace getfem {
Normal, m);
GMM_ASSERT1(local_workspace.assembled_tensor().size() == m.dim(),
- "Wrong dimension of the tranformation expression");
+ "Wrong dimension of the transformation expression");
P.resize(m.dim());
gmm::copy(local_workspace.assembled_tensor().as_vector(), P);
diff --git a/src/getfem_level_set.cc b/src/getfem_level_set.cc
index b30bea8..2b17804 100644
--- a/src/getfem_level_set.cc
+++ b/src/getfem_level_set.cc
@@ -25,7 +25,7 @@
namespace getfem {
level_set::level_set(const mesh &msh, dim_type deg,
- bool with_secondary_)
+ bool with_secondary_)
: degree_(deg), mf(&classical_mesh_fem(msh, deg)),
with_secondary(with_secondary_) {
shift_ls = scalar_type(0);
@@ -64,20 +64,20 @@ namespace getfem {
}
pmesher_signed_distance level_set::mls_of_convex(size_type cv, unsigned
lsnum,
- bool inverted) const {
+ bool inverted) const {
assert(this); assert(mf);
GMM_ASSERT1(mf->linked_mesh().convex_index().is_in(cv), "convex " << cv
- << " is not in the level set mesh!");
+ << " is not in the level set mesh!");
GMM_ASSERT1(mf->fem_of_element(cv), "Internal error");
GMM_ASSERT1(!mf->is_reduced(), "Internal error");
std::vector<scalar_type> coeff(mf->nb_basic_dof_of_element(cv));
GMM_ASSERT1(values(lsnum).size() == mf->nb_dof(),
- "Inconsistent state in the levelset: nb_dof=" <<
- mf->nb_dof() << ", values(" << lsnum << ").size=" <<
- values(lsnum).size());
+ "Inconsistent state in the levelset: nb_dof=" <<
+ mf->nb_dof() << ", values(" << lsnum << ").size=" <<
+ values(lsnum).size());
for (size_type i = 0; i < coeff.size(); ++i)
coeff[i] = (!inverted ? scalar_type(1) : scalar_type(-1)) *
- values(lsnum)[mf->ind_basic_dof_of_element(cv)[i]];
+ values(lsnum)[mf->ind_basic_dof_of_element(cv)[i]];
return new_mesher_level_set(mf->fem_of_element(cv), coeff, shift_ls);
}
@@ -89,16 +89,16 @@ namespace getfem {
void level_set::simplify(scalar_type eps) {
for (dal::bv_visitor cv(mf->linked_mesh().convex_index());
- !cv.finished(); ++cv) {
+ !cv.finished(); ++cv) {
scalar_type h = mf->linked_mesh().convex_radius_estimate(cv);
for (size_type i = 0; i < mf->nb_basic_dof_of_element(cv); ++i) {
- size_type dof = mf->ind_basic_dof_of_element(cv)[i];
- if (gmm::abs(primary_[dof]) < h*eps) {
- primary_[dof] = scalar_type(0);
- // cout << "Simplify dof " << dof << " : " << mf->point_of_dof(dof)
<< endl;
- }
- if (has_secondary() && gmm::abs(secondary_[dof]) < h*eps)
- secondary_[dof] = scalar_type(0);
+ size_type dof = mf->ind_basic_dof_of_element(cv)[i];
+ if (gmm::abs(primary_[dof]) < h*eps) {
+ primary_[dof] = scalar_type(0);
+ // cout << "Simplify dof " << dof << " : " << mf->point_of_dof(dof)
<< endl;
+ }
+ if (has_secondary() && gmm::abs(secondary_[dof]) < h*eps)
+ secondary_[dof] = scalar_type(0);
}
}
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index 69aae4b..ee591f6 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -66,7 +66,7 @@ namespace getfem {
if (r[ic][0])
for (size_type jc = 0; jc < icv.size(); ++jc)
r.add(icv[jc]);
-
+
bgeot::geotrans_inv_convex giv;
for (short_type f = 0; f < pgt->structure()->nb_faces(); ++f) {
if (r[ic][f+1]) {
@@ -86,17 +86,17 @@ namespace getfem {
base_node pt, barycentre
=gmm::mean_value(pgtsub->convex_ref()->points_of_face(fsub));
pt = pgtsub->transform(barycentre, points_of_convex(icv[jc]));
-
+
giv.init(points_of_convex(ic), pgt);
giv.invert(pt, barycentre);
-
+
if (gmm::abs(pgt->convex_ref()->is_in_face(f,barycentre)) <
0.001)
r.add(icv[jc], fsub);
}
}
}
- }
+ }
}
for (size_type jc = 0; jc < icv.size(); ++jc)
@@ -111,10 +111,10 @@ namespace getfem {
base_node pt, barycentre
= gmm::mean_value(pgtsub->convex_ref()->points_of_face(fsub));
pt = pgtsub->transform(barycentre, points_of_convex(icv[jc]));
-
+
giv.init(points_of_convex(ic), pgt);
giv.invert(pt, barycentre);
-
+
for (short_type f = 0; f < pgt->structure()->nb_faces(); ++f)
if (gmm::abs(pgt->convex_ref()->is_in_face(f,barycentre)) <
0.001)
{ r.add(ic, f); break; }
@@ -136,7 +136,7 @@ namespace getfem {
: bgeot::basic_mesh(m), name_(name) { init(); }
mesh::mesh(const mesh &m) : context_dependencies(),
- std::enable_shared_from_this<getfem::mesh>()
+ std::enable_shared_from_this<getfem::mesh>()
{ copy_from(m); }
mesh &mesh::operator=(const mesh &m) {
@@ -148,58 +148,58 @@ namespace getfem {
#if GETFEM_PARA_LEVEL > 1
- void mesh::compute_mpi_region(void) const {
- int size, rank;
- MPI_Comm_rank(MPI_COMM_WORLD, &rank);
- MPI_Comm_size(MPI_COMM_WORLD, &size);
-
- if (size < 2) {
- mpi_region = mesh_region::all_convexes();
- mpi_region.from_mesh(*this);
- } else {
- int ne = int(nb_convex());
- std::vector<int> xadj(ne+1), adjncy, numelt(ne), npart(ne);
- std::vector<int> indelt(nb_allocated_convex());
-
- double t_ref = MPI_Wtime();
-
- int j = 0, k = 0;
- ind_set s;
- for (dal::bv_visitor ic(convex_index()); !ic.finished(); ++ic, ++j) {
- numelt[j] = ic;
- indelt[ic] = j;
- }
- j = 0;
- for (dal::bv_visitor ic(convex_index()); !ic.finished(); ++ic, ++j) {
- xadj[j] = k;
- neighbours_of_convex(ic, s);
- for (ind_set::iterator it = s.begin();
- it != s.end(); ++it) { adjncy.push_back(indelt[*it]); ++k; }
- }
+ void mesh::compute_mpi_region() const {
+ int size, rank;
+ MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+ MPI_Comm_size(MPI_COMM_WORLD, &size);
+
+ if (size < 2) {
+ mpi_region = mesh_region::all_convexes();
+ mpi_region.from_mesh(*this);
+ } else {
+ int ne = int(nb_convex());
+ std::vector<int> xadj(ne+1), adjncy, numelt(ne), npart(ne);
+ std::vector<int> indelt(nb_allocated_convex());
+
+ double t_ref = MPI_Wtime();
+
+ int j = 0, k = 0;
+ ind_set s;
+ for (dal::bv_visitor ic(convex_index()); !ic.finished(); ++ic, ++j) {
+ numelt[j] = ic;
+ indelt[ic] = j;
+ }
+ j = 0;
+ for (dal::bv_visitor ic(convex_index()); !ic.finished(); ++ic, ++j) {
xadj[j] = k;
+ neighbours_of_convex(ic, s);
+ for (ind_set::iterator it = s.begin();
+ it != s.end(); ++it) { adjncy.push_back(indelt[*it]); ++k; }
+ }
+ xadj[j] = k;
#ifdef GETFEM_HAVE_METIS_OLD_API
- int wgtflag = 0, numflag = 0, edgecut;
- int options[5] = {0,0,0,0,0};
- METIS_PartGraphKway(&ne, &(xadj[0]), &(adjncy[0]), 0, 0, &wgtflag,
- &numflag, &size, options, &edgecut, &(npart[0]));
+ int wgtflag = 0, numflag = 0, edgecut;
+ int options[5] = {0,0,0,0,0};
+ METIS_PartGraphKway(&ne, &(xadj[0]), &(adjncy[0]), 0, 0, &wgtflag,
+ &numflag, &size, options, &edgecut, &(npart[0]));
#else
- int ncon = 1, edgecut;
- int options[METIS_NOPTIONS] = { 0 };
- METIS_SetDefaultOptions(options);
- METIS_PartGraphKway(&ne, &ncon, &(xadj[0]), &(adjncy[0]), 0, 0, 0, &size,
- 0, 0, options, &edgecut, &(npart[0]));
+ int ncon = 1, edgecut;
+ int options[METIS_NOPTIONS] = { 0 };
+ METIS_SetDefaultOptions(options);
+ METIS_PartGraphKway(&ne, &ncon, &(xadj[0]), &(adjncy[0]), 0, 0, 0,
+ &size, 0, 0, options, &edgecut, &(npart[0]));
#endif
- for (size_type i = 0; i < size_type(ne); ++i)
- if (npart[i] == rank) mpi_region.add(numelt[i]);
+ for (size_type i = 0; i < size_type(ne); ++i)
+ if (npart[i] == rank) mpi_region.add(numelt[i]);
- if (MPI_IS_MASTER())
- cout << "Partition time "<< MPI_Wtime()-t_ref << endl;
- }
- modified = false;
- valid_sub_regions.clear();
+ if (MPI_IS_MASTER())
+ cout << "Partition time "<< MPI_Wtime()-t_ref << endl;
}
+ modified = false;
+ valid_sub_regions.clear();
+ }
void mesh::compute_mpi_sub_region(size_type n) const {
if (valid_cvf_sets.is_in(n)) {
@@ -236,15 +236,15 @@ namespace getfem {
if (with_renumbering) { // Could be optimized no using only swap_convex
std::vector<size_type> cmk, iord(nbc), iordinv(nbc);
for (i = 0; i < nbc; ++i) iord[i] = iordinv[i] = i;
-
+
bgeot::cuthill_mckee_on_convexes(*this, cmk);
for (i = 0; i < nbc; ++i) {
- j = iordinv[cmk[i]];
- if (i != j) {
- swap_convex(i, j);
- std::swap(iord[i], iord[j]);
- std::swap(iordinv[iord[i]], iordinv[iord[j]]);
- }
+ j = iordinv[cmk[i]];
+ if (i != j) {
+ swap_convex(i, j);
+ std::swap(iord[i], iord[j]);
+ std::swap(iordinv[iord[i]], iordinv[iord[j]]);
+ }
}
}
}
@@ -433,7 +433,7 @@ namespace getfem {
}
return r;
}
-
+
scalar_type mesh::maximal_convex_radius_estimate() const {
if (convex_index().empty()) return 1;
scalar_type r = convex_radius_estimate(convex_index().first_true());
@@ -487,29 +487,29 @@ namespace getfem {
{ te = true; }
else if (!bgeot::casecmp(tmp, "POINT")) {
bgeot::get_token(ist, tmp);
- if (!bgeot::casecmp(tmp, "COUNT")) {
- bgeot::get_token(ist, tmp); // Ignored. Used in some applications
- } else {
- size_type ip = atoi(tmp.c_str());
- dim_type d = 0;
- GMM_ASSERT1(!npt.is_in(ip),
- "Two points with the same index. loading aborted.");
- npt.add(ip);
- bgeot::get_token(ist, tmp);
- while (isdigit(tmp[0]) || tmp[0] == '-' || tmp[0] == '+'
- || tmp[0] == '.')
- { tmpv[d++] = atof(tmp.c_str()); bgeot::get_token(ist, tmp); }
- please_get = false;
- base_node v(d);
- for (size_type i = 0; i < d; i++) v[i] = tmpv[i];
- size_type ipl = add_point(v);
- if (ip != ipl) {
- GMM_ASSERT1(!npt.is_in(ipl), "Two points [#" << ip << " and #"
- << ipl << "] with the same coords "<< v
- << ". loading aborted.");
- swap_points(ip, ipl);
- }
- }
+ if (!bgeot::casecmp(tmp, "COUNT")) {
+ bgeot::get_token(ist, tmp); // Ignored. Used in some applications
+ } else {
+ size_type ip = atoi(tmp.c_str());
+ dim_type d = 0;
+ GMM_ASSERT1(!npt.is_in(ip),
+ "Two points with the same index. loading aborted.");
+ npt.add(ip);
+ bgeot::get_token(ist, tmp);
+ while (isdigit(tmp[0]) || tmp[0] == '-' || tmp[0] == '+'
+ || tmp[0] == '.')
+ { tmpv[d++] = atof(tmp.c_str()); bgeot::get_token(ist, tmp); }
+ please_get = false;
+ base_node v(d);
+ for (size_type i = 0; i < d; i++) v[i] = tmpv[i];
+ size_type ipl = add_point(v);
+ if (ip != ipl) {
+ GMM_ASSERT1(!npt.is_in(ipl), "Two points [#" << ip << " and #"
+ << ipl << "] with the same coords "<< v
+ << ". loading aborted.");
+ swap_points(ip, ipl);
+ }
+ }
} else if (tmp.size()) {
GMM_ASSERT1(false, "Syntax error in file, at token '" << tmp
<< "', pos=" << std::streamoff(ist.tellg()));
@@ -533,30 +533,30 @@ namespace getfem {
else if (!bgeot::casecmp(tmp, "CONVEX")) {
size_type ic;
bgeot::get_token(ist, tmp);
- if (!bgeot::casecmp(tmp, "COUNT")) {
- bgeot::get_token(ist, tmp); // Ignored. Used in some applications
- } else {
- ic = gmm::abs(atoi(tmp.c_str()));
- GMM_ASSERT1(!ncv.is_in(ic),
- "Negative or repeated index, loading aborted.");
- ncv.add(ic);
-
- int rgt = bgeot::get_token(ist, tmp);
- if (rgt != 3) { // for backward compatibility with version 1.7
- char c; ist.get(c);
- while (!isspace(c)) { tmp.push_back(c); ist.get(c); }
- }
-
- bgeot::pgeometric_trans pgt = bgeot::geometric_trans_descriptor(tmp);
- size_type nb = pgt->nb_points();
-
- cv[ic].cstruct = pgt;
- cv[ic].pts.resize(nb);
- for (size_type i = 0; i < nb; i++) {
- bgeot::get_token(ist, tmp);
- cv[ic].pts[i] = gmm::abs(atoi(tmp.c_str()));
- }
- }
+ if (!bgeot::casecmp(tmp, "COUNT")) {
+ bgeot::get_token(ist, tmp); // Ignored. Used in some applications
+ } else {
+ ic = gmm::abs(atoi(tmp.c_str()));
+ GMM_ASSERT1(!ncv.is_in(ic),
+ "Negative or repeated index, loading aborted.");
+ ncv.add(ic);
+
+ int rgt = bgeot::get_token(ist, tmp);
+ if (rgt != 3) { // for backward compatibility with version 1.7
+ char c; ist.get(c);
+ while (!isspace(c)) { tmp.push_back(c); ist.get(c); }
+ }
+
+ bgeot::pgeometric_trans pgt = bgeot::geometric_trans_descriptor(tmp);
+ size_type nb = pgt->nb_points();
+
+ cv[ic].cstruct = pgt;
+ cv[ic].pts.resize(nb);
+ for (size_type i = 0; i < nb; i++) {
+ bgeot::get_token(ist, tmp);
+ cv[ic].pts[i] = gmm::abs(atoi(tmp.c_str()));
+ }
+ }
}
else if (tmp.size()) {
GMM_ASSERT1(false, "Syntax error reading a mesh file "
@@ -748,7 +748,7 @@ namespace getfem {
a transformation of any pgt to the equivalent equilateral pgt
(for prisms etc) */
if (bgeot::basic_structure(pgt->structure())
- == bgeot::simplex_structure(P))
+ == bgeot::simplex_structure(P))
gmm::mult(base_matrix(K),equilateral_to_GT_PK_grad(P),K);
q = std::max(q, gmm::condition_number(K));
}
@@ -823,23 +823,23 @@ namespace getfem {
mr.error_if_not_convexes();
dal::bit_vector visited;
bgeot::mesh_structure::ind_set neighbours;
-
+
for (mr_visitor i(mr); !i.finished(); ++i) {
size_type cv1 = i.cv();
short_type nbf = m.structure_of_convex(i.cv())->nb_faces();
bool neighbour_visited = false;
for (short_type f = 0; f < nbf; ++f) {
- neighbours.resize(0); m.neighbours_of_convex(cv1, f, neighbours);
- for (size_type j = 0; j < neighbours.size(); ++j)
- if (visited.is_in(neighbours[j]))
- { neighbour_visited = true; break; }
+ neighbours.resize(0); m.neighbours_of_convex(cv1, f, neighbours);
+ for (size_type j = 0; j < neighbours.size(); ++j)
+ if (visited.is_in(neighbours[j]))
+ { neighbour_visited = true; break; }
}
if (!neighbour_visited) {
for (short_type f = 0; f < nbf; ++f) {
- size_type cv2 = m.neighbour_of_convex(cv1, f);
- if (cv2 != size_type(-1) && mr.is_in(cv2)) mrr.add(cv1,f);
- }
- visited.add(cv1);
+ size_type cv2 = m.neighbour_of_convex(cv1, f);
+ if (cv2 != size_type(-1) && mr.is_in(cv2)) mrr.add(cv1,f);
+ }
+ visited.add(cv1);
}
}
@@ -847,16 +847,16 @@ namespace getfem {
size_type cv1 = i.cv();
if (!(visited.is_in(cv1))) {
short_type nbf = m.structure_of_convex(i.cv())->nb_faces();
- for (short_type f = 0; f < nbf; ++f) {
- neighbours.resize(0); m.neighbours_of_convex(cv1, f, neighbours);
- bool ok = false;
- for (size_type j = 0; j < neighbours.size(); ++j) {
- if (visited.is_in(neighbours[j])) { ok = false; break; }
- if (mr.is_in(neighbours[j])) { ok = true; }
- }
- if (ok) { mrr.add(cv1,f); }
- }
- visited.add(cv1);
+ for (short_type f = 0; f < nbf; ++f) {
+ neighbours.resize(0); m.neighbours_of_convex(cv1, f, neighbours);
+ bool ok = false;
+ for (size_type j = 0; j < neighbours.size(); ++j) {
+ if (visited.is_in(neighbours[j])) { ok = false; break; }
+ if (mr.is_in(neighbours[j])) { ok = true; }
+ }
+ if (ok) { mrr.add(cv1,f); }
+ }
+ visited.add(cv1);
}
}
return mrr;
@@ -881,7 +881,7 @@ namespace getfem {
const base_node &pt1, const base_node &pt2) {
mesh_region mrr;
size_type N = m.dim();
- GMM_ASSERT1(pt1.size() == N && pt2.size() == N, "Wrong dimensions");
+ GMM_ASSERT1(pt1.size() == N && pt2.size() == N, "Wrong dimensions");
for (getfem::mr_visitor i(mr, m); !i.finished(); ++i)
if (i.f() != short_type(-1)) {
bgeot::mesh_structure::ind_pt_face_ct pt
@@ -904,7 +904,7 @@ namespace getfem {
const base_node &pt1, const base_node
&pt2) {
mesh_region mrr;
size_type N = m.dim();
- GMM_ASSERT1(pt1.size() == N && pt2.size() == N, "Wrong dimensions");
+ GMM_ASSERT1(pt1.size() == N && pt2.size() == N, "Wrong dimensions");
for (getfem::mr_visitor i(mr, m); !i.finished(); ++i)
if (i.f() == short_type(-1)) {
bgeot::mesh_structure::ind_cv_ct pt = m.ind_points_of_convex(i.cv());
diff --git a/src/getfem_mesh_im_level_set.cc b/src/getfem_mesh_im_level_set.cc
index 6fd7afd..9b0bef4 100644
--- a/src/getfem_mesh_im_level_set.cc
+++ b/src/getfem_mesh_im_level_set.cc
@@ -742,8 +742,6 @@ namespace getfem {
is_adapted = true; touch();
}
-
-
} /* end of namespace getfem. */
diff --git a/src/getfem_mesh_region.cc b/src/getfem_mesh_region.cc
index 11ebfdc..3422010 100644
--- a/src/getfem_mesh_region.cc
+++ b/src/getfem_mesh_region.cc
@@ -35,9 +35,9 @@ namespace getfem {
mesh_region::mesh_region() : p(std::make_shared<impl>()), id_(size_type(-2)),
- type_(size_type(-1)),
+ type_(size_type(-1)),
partitioning_allowed(true), parent_mesh(0), index_updated(false)
- {
+ {
if (me_is_multithreaded_now()) prohibit_partitioning();
}
@@ -45,36 +45,36 @@ namespace getfem {
partitioning_allowed(true), parent_mesh(0), index_updated(false)
{ }
- mesh_region::mesh_region(mesh& m, size_type id__, size_type type) :
+ mesh_region::mesh_region(mesh& m, size_type id__, size_type type) :
p(std::make_shared<impl>()), id_(id__), type_(type),
partitioning_allowed(true), parent_mesh(&m),
index_updated(false)
- {
- if (me_is_multithreaded_now()) prohibit_partitioning();
+ {
+ if (me_is_multithreaded_now()) prohibit_partitioning();
}
- mesh_region::mesh_region(const dal::bit_vector &bv) :
+ mesh_region::mesh_region(const dal::bit_vector &bv) :
p(std::make_shared<impl>()), id_(size_type(-2)), type_(size_type(-1)),
partitioning_allowed(true), parent_mesh(0), index_updated(false)
- {
- if (me_is_multithreaded_now()) prohibit_partitioning();
- add(bv);
+ {
+ if (me_is_multithreaded_now()) prohibit_partitioning();
+ add(bv);
}
- void mesh_region::touch_parent_mesh()
+ void mesh_region::touch_parent_mesh()
{
if (parent_mesh) parent_mesh->touch_from_region(id_);
}
- const mesh_region& mesh_region::from_mesh(const mesh &m) const
+ const mesh_region& mesh_region::from_mesh(const mesh &m) const
{
- if (!p.get())
+ if (!p.get())
{
mesh_region *r = const_cast<mesh_region*>(this);
- if (id_ == size_type(-1))
+ if (id_ == size_type(-1))
{
r->p = std::make_shared<impl>();
r->add(m.convex_index());
- }
+ }
else if (id_ != size_type(-2))
{
*r = m.region(id_);
@@ -84,7 +84,7 @@ namespace getfem {
return *this;
}
- mesh_region& mesh_region::operator=(const mesh_region &from)
+ mesh_region& mesh_region::operator=(const mesh_region &from)
{
if (!this->parent_mesh && !from.parent_mesh)
{
@@ -137,7 +137,7 @@ namespace getfem {
return true;
}
- face_bitset mesh_region::operator[](size_t cv) const
+ face_bitset mesh_region::operator[](size_t cv) const
{
map_t::const_iterator it = rp().m.find(cv);
if (it != rp().m.end()) return (*it).second;
@@ -151,11 +151,11 @@ namespace getfem {
itend = partition_end ();
index_updated = true;
}
- mesh_region::const_iterator
+ mesh_region::const_iterator
mesh_region::partition_begin( ) const
{
size_type region_size = rp().m.size();
- if (region_size < num_threads())
+ if (region_size < num_threads())
{ //for small regions: put the whole region into zero thread
if (this_thread() == 0) return rp().m.begin();
else return rp().m.end();
@@ -171,11 +171,11 @@ namespace getfem {
return it;
}
- mesh_region::const_iterator
+ mesh_region::const_iterator
mesh_region::partition_end( ) const
{
size_type region_size = rp().m.size();
- if (region_size < num_threads()) return rp().m.end();
+ if (region_size < num_threads()) return rp().m.end();
size_type partition_size = static_cast<size_type>
(std::ceil(static_cast<scalar_type>(region_size)/
@@ -188,21 +188,21 @@ namespace getfem {
return it;
}
- mesh_region::const_iterator mesh_region::begin( ) const
+ mesh_region::const_iterator mesh_region::begin( ) const
{
GMM_ASSERT1(p != 0, "Internal error");
- if (me_is_multithreaded_now() && partitioning_allowed)
- {
+ if (me_is_multithreaded_now() && partitioning_allowed)
+ {
update_partition_iterators();
return itbegin;
}
else { return rp().m.begin(); }
}
- mesh_region::const_iterator mesh_region::end ( ) const
- {
- if (me_is_multithreaded_now() && partitioning_allowed)
- {
+ mesh_region::const_iterator mesh_region::end ( ) const
+ {
+ if (me_is_multithreaded_now() && partitioning_allowed)
+ {
update_partition_iterators();
return itend;
}
@@ -239,18 +239,18 @@ namespace getfem {
}
/* may be optimized .. */
- const dal::bit_vector& mesh_region::index() const
+ const dal::bit_vector& mesh_region::index() const
{
dal::bit_vector& convex_index = rp().index_.thrd_cast();
convex_index.clear();
- for (const_iterator it = begin(); it != end(); ++it)
+ for (const_iterator it = begin(); it != end(); ++it)
{
if (it->second.any()) convex_index.add(it->first);
}
return convex_index;
}
- void mesh_region::add(const dal::bit_vector &bv)
+ void mesh_region::add(const dal::bit_vector &bv)
{
for (dal::bv_visitor i(bv); !i.finished(); ++i)
{
@@ -260,15 +260,15 @@ namespace getfem {
index_updated = false;
}
- void mesh_region::add(size_type cv, short_type f)
+ void mesh_region::add(size_type cv, short_type f)
{
- wp().m[cv].set(short_type(f+1),1);
+ wp().m[cv].set(short_type(f+1),1);
touch_parent_mesh();
index_updated = false;
}
- void mesh_region::sup_all(size_type cv)
- {
+ void mesh_region::sup_all(size_type cv)
+ {
map_t::iterator it = wp().m.find(cv);
if (it != wp().m.end()) {
wp().m.erase(it);
@@ -277,29 +277,29 @@ namespace getfem {
index_updated = false;
}
- void mesh_region::sup(size_type cv, short_type f)
- {
+ void mesh_region::sup(size_type cv, short_type f)
+ {
map_t::iterator it = wp().m.find(cv);
if (it != wp().m.end()) {
(*it).second.set(short_type(f+1),0);
- if ((*it).second.none()) wp().m.erase(it);
+ if ((*it).second.none()) wp().m.erase(it);
touch_parent_mesh();
}
index_updated = false;
}
- void mesh_region::clear()
- {
+ void mesh_region::clear()
+ {
wp().m.clear(); touch_parent_mesh();
index_updated = false;
}
- void mesh_region::clean()
+ void mesh_region::clean()
{
- for (map_t::iterator it = wp().m.begin(), itn;
- it != wp().m.end(); it = itn)
+ for (map_t::iterator it = wp().m.begin(), itn;
+ it != wp().m.end(); it = itn)
{
- itn = it;
+ itn = it;
++itn;
if ( !(*it).second.any() ) wp().m.erase(it);
}
@@ -324,7 +324,7 @@ namespace getfem {
index_updated = false;
}
- bool mesh_region::is_in(size_type cv, short_type f) const
+ bool mesh_region::is_in(size_type cv, short_type f) const
{
GMM_ASSERT1(p.get(), "Use from mesh on that region before");
map_t::const_iterator it = rp().m.find(cv);
@@ -332,12 +332,12 @@ namespace getfem {
return ((*it).second)[short_type(f+1)];
}
- bool mesh_region::is_in(size_type cv, short_type f, const mesh &m) const
+ bool mesh_region::is_in(size_type cv, short_type f, const mesh &m) const
{
if (p.get()) {
map_t::const_iterator it = rp().m.find(cv);
if (it == rp().m.end() || short_type(f+1) >= MAX_FACES_PER_CV)
- return false;
+ return false;
return ((*it).second)[short_type(f+1)];
}
else
@@ -347,34 +347,34 @@ namespace getfem {
}
}
-
- bool mesh_region::is_empty() const
+
+ bool mesh_region::is_empty() const
{
return rp().m.empty();
}
- bool mesh_region::is_only_convexes() const
+ bool mesh_region::is_only_convexes() const
{
- return is_empty() ||
- (or_mask()[0] == true && or_mask().count() == 1);
+ return is_empty() ||
+ (or_mask()[0] == true && or_mask().count() == 1);
}
- bool mesh_region::is_only_faces() const
- {
- return is_empty() || (and_mask()[0] == false);
+ bool mesh_region::is_only_faces() const
+ {
+ return is_empty() || (and_mask()[0] == false);
}
- face_bitset mesh_region::faces_of_convex(size_type cv) const
+ face_bitset mesh_region::faces_of_convex(size_type cv) const
{
map_t::const_iterator it = rp().m.find(cv);
- if (it != rp().m.end()) return ((*it).second) >> 1;
+ if (it != rp().m.end()) return ((*it).second) >> 1;
else return face_bitset();
}
- face_bitset mesh_region::and_mask() const
+ face_bitset mesh_region::and_mask() const
{
- face_bitset bs;
+ face_bitset bs;
if (rp().m.empty()) return bs;
bs.set();
for (map_t::const_iterator it = rp().m.begin(); it != rp().m.end(); ++it)
@@ -382,16 +382,16 @@ namespace getfem {
return bs;
}
- face_bitset mesh_region::or_mask() const
+ face_bitset mesh_region::or_mask() const
{
- face_bitset bs;
+ face_bitset bs;
if (rp().m.empty()) return bs;
for (map_t::const_iterator it = rp().m.begin(); it != rp().m.end(); ++it)
if ( (*it).second.any() ) bs |= (*it).second;
return bs;
}
- size_type mesh_region::size() const
+ size_type mesh_region::size() const
{
size_type sz=0;
for (map_t::const_iterator it = begin(); it != end(); ++it)
@@ -399,7 +399,7 @@ namespace getfem {
return sz;
}
- size_type mesh_region::unpartitioned_size() const
+ size_type mesh_region::unpartitioned_size() const
{
size_type sz=0;
for (map_t::const_iterator it = rp().m.begin(); it != rp().m.end(); ++it)
@@ -408,8 +408,8 @@ namespace getfem {
}
- mesh_region mesh_region::intersection(const mesh_region &a,
- const mesh_region &b)
+ mesh_region mesh_region::intersection(const mesh_region &a,
+ const mesh_region &b)
{
GMM_TRACE4("intersection of "<<a.id()<<" and "<<b.id());
mesh_region r;
@@ -418,20 +418,20 @@ namespace getfem {
(they only exist to provide a default argument to the mesh_region
parameters of assembly procedures etc. */
GMM_ASSERT1(a.id() != size_type(-1)||
- b.id() != size_type(-1), "the 'all_convexes' regions "
- "are not supported for set operations");
- if (a.id() == size_type(-1))
+ b.id() != size_type(-1), "the 'all_convexes' regions "
+ "are not supported for set operations");
+ if (a.id() == size_type(-1))
{
for (const_iterator it = b.begin();it != b.end(); ++it)
r.wp().m.insert(*it);
return r;
}
- else if (b.id() == size_type(-1))
+ else if (b.id() == size_type(-1))
{
for (const_iterator it = a.begin();it != a.end(); ++it)
r.wp().m.insert(*it);
return r;
}
- map_t::const_iterator
+ map_t::const_iterator
ita = a.begin(), enda = a.end(),
itb = b.begin(), endb = b.end();
@@ -450,8 +450,8 @@ namespace getfem {
return r;
}
- mesh_region mesh_region::merge(const mesh_region &a,
- const mesh_region &b)
+ mesh_region mesh_region::merge(const mesh_region &a,
+ const mesh_region &b)
{
GMM_TRACE4("Merger of " << a.id() << " and " << b.id());
mesh_region r;
@@ -459,7 +459,7 @@ namespace getfem {
b.id() != size_type(-1), "the 'all_convexes' regions "
"are not supported for set operations");
for (const_iterator it = a.begin();it != a.end(); ++it)
- {
+ {
r.wp().m.insert(*it);
}
for (const_iterator it = b.begin();it != b.end(); ++it)
@@ -471,38 +471,38 @@ namespace getfem {
mesh_region mesh_region::subtract(const mesh_region &a,
- const mesh_region &b)
+ const mesh_region &b)
{
GMM_TRACE4("subtraction of "<<a.id()<<" and "<<b.id());
mesh_region r;
GMM_ASSERT1(a.id() != size_type(-1) &&
b.id() != size_type(-1), "the 'all_convexes' regions "
"are not supported for set operations");
- for (const_iterator ita = a.begin();ita != a.end(); ++ita)
+ for (const_iterator ita = a.begin();ita != a.end(); ++ita)
r.wp().m.insert(*ita);
- for (const_iterator itb = b.begin();itb != b.end(); ++itb)
+ for (const_iterator itb = b.begin();itb != b.end(); ++itb)
{
size_type cv = itb->first;
map_t::iterator it = r.wp().m.find(cv);
if (it != r.wp().m.end()){
- it->second &= ~(itb->second);
- if (it->second.none()) r.wp().m.erase(it);
- }
+ it->second &= ~(itb->second);
+ if (it->second.none()) r.wp().m.erase(it);
+ }
}
return r;
}
int mesh_region::region_is_faces_of(const getfem::mesh& m1,
- const mesh_region &rg2,
- const getfem::mesh& m2) const {
+ const mesh_region &rg2,
+ const getfem::mesh& m2) const {
if (&m1 != &m2) return 0;
int r = 1, partially = 0;
for (mr_visitor cv(*this, m1); !cv.finished(); cv.next())
if (cv.is_face() && rg2.is_in(cv.cv(),short_type(-1), m2))
- partially = -1;
+ partially = -1;
else
- r = 0;
+ r = 0;
if (r == 1) return 1; else return partially;
}
@@ -512,17 +512,17 @@ namespace getfem {
}
- void mesh_region::error_if_not_faces() const
+ void mesh_region::error_if_not_faces() const
{
GMM_ASSERT1(is_only_faces(), "Expecting a set of faces, not convexes");
}
- void mesh_region::error_if_not_convexes() const
+ void mesh_region::error_if_not_convexes() const
{
GMM_ASSERT1(is_only_convexes(), "Expecting a set of convexes, not faces");
}
- void mesh_region::error_if_not_homogeneous() const
+ void mesh_region::error_if_not_homogeneous() const
{
GMM_ASSERT1(is_only_faces() || is_only_convexes(), "Expecting a set "
"of convexes or a set of faces, but not a mixed set");
@@ -532,60 +532,60 @@ namespace getfem {
#if GETFEM_PARA_LEVEL > 1
-
+
mesh_region::visitor::visitor(const mesh_region &s, const mesh &m,
- bool intersect_with_mpi) :
- cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
+ bool intersect_with_mpi) :
+ cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
{
if ((me_is_multithreaded_now() && s.partitioning_allowed)) {
s.from_mesh(m);
init(s);
} else {
if (s.id() == size_type(-1)) {
- if (intersect_with_mpi)
- init(m.get_mpi_region());
- else
- init(m.convex_index());
+ if (intersect_with_mpi)
+ init(m.get_mpi_region());
+ else
+ init(m.convex_index());
} else if (s.p.get()) {
- if (intersect_with_mpi) {
- mpi_rg = std::make_unique<mesh_region>(s);
- mpi_rg->from_mesh(m);
- m.intersect_with_mpi_region(*mpi_rg);
- init(*mpi_rg);
- } else
- init(s);
+ if (intersect_with_mpi) {
+ mpi_rg = std::make_unique<mesh_region>(s);
+ mpi_rg->from_mesh(m);
+ m.intersect_with_mpi_region(*mpi_rg);
+ init(*mpi_rg);
+ } else
+ init(s);
} else {
- if (intersect_with_mpi)
- init(m.get_mpi_sub_region(s.id()));
- else
- init(m.region(s.id()));
+ if (intersect_with_mpi)
+ init(m.get_mpi_sub_region(s.id()));
+ else
+ init(m.region(s.id()));
}
}
}
#else
-
+
mesh_region::visitor::visitor(const mesh_region &s, const mesh &m, bool)
- :cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
+ :cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
{
if ((me_is_multithreaded_now() && s.partitioning_allowed)) {
s.from_mesh(m);
init(s);
} else {
if (s.id() == size_type(-1)) {
- init(m.convex_index());
+ init(m.convex_index());
} else if (s.p.get()) {
- init(s);
+ init(s);
} else {
- init(m.region(s.id()));
+ init(m.region(s.id()));
}
}
}
#endif
-
- bool mesh_region::visitor::next()
+
+ bool mesh_region::visitor::next()
{
if (whole_mesh) {
if (itb == iteb) { finished_ = true; return false; }
@@ -595,26 +595,26 @@ namespace getfem {
++itb; while (itb != iteb && !(*itb)) ++itb;
return true;
}
- while (c.none())
+ while (c.none())
{
- if (it == ite) { finished_=true; return false; }
- cv_ = it->first;
- c = it->second;
- f_ = short_type(-1);
- ++it;
- if (c.none()) continue;
+ if (it == ite) { finished_=true; return false; }
+ cv_ = it->first;
+ c = it->second;
+ f_ = short_type(-1);
+ ++it;
+ if (c.none()) continue;
}
next_face();
return true;
}
mesh_region::visitor::visitor(const mesh_region &s) :
- cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
+ cv_(size_type(-1)), f_(short_type(-1)), finished_(false)
{
init(s);
}
- void mesh_region::visitor::init(const dal::bit_vector &bv)
+ void mesh_region::visitor::init(const dal::bit_vector &bv)
{
whole_mesh = true;
itb = bv.begin(); iteb = bv.end();
@@ -622,7 +622,7 @@ namespace getfem {
next();
}
- void mesh_region::visitor::init(const mesh_region &s)
+ void mesh_region::visitor::init(const mesh_region &s)
{
whole_mesh = false;
it = s.begin();
@@ -630,18 +630,18 @@ namespace getfem {
next();
}
- std::ostream & operator <<(std::ostream &os, const mesh_region &w)
+ std::ostream & operator <<(std::ostream &os, const mesh_region &w)
{
if (w.id() == size_type(-1))
os << " ALL_CONVEXES";
else if (w.p.get())
{
- for (mr_visitor cv(w); !cv.finished(); cv.next())
- {
- os << cv.cv();
- if (cv.is_face()) os << "/" << cv.f();
- os << " ";
- }
+ for (mr_visitor cv(w); !cv.finished(); cv.next())
+ {
+ os << cv.cv();
+ if (cv.is_face()) os << "/" << cv.f();
+ os << " ";
+ }
}
else
{
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index c2c7e81..f7eaf97 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -3917,9 +3917,8 @@ model_complex_plain_vector &
size_type ib = add_linear_term
(md, mim, expr, region, true, true, "Generic elliptic", true);
if (ib == size_type(-1))
- ib = add_nonlinear_term
- (md, mim, expr, region, false, false,
- "Generic elliptic (nonlinear)");
+ ib = add_nonlinear_term(md, mim, expr, region, false, false,
+ "Generic elliptic (nonlinear)");
return ib;
}
}
@@ -4075,9 +4074,8 @@ model_complex_plain_vector &
size_type ib = add_source_term_generic_assembly_brick
(md, mim, expr, region, "Source term", varname, directdataname, true);
if (ib == size_type(-1)) {
- ib = add_nonlinear_term
- (md, mim, "-("+expr+")", region, false, false,
- "Source term (nonlinear)");
+ ib = add_nonlinear_term(md, mim, "-("+expr+")", region, false, false,
+ "Source term (nonlinear)");
if (directdataname.size())
add_source_term_generic_assembly_brick
(md, mim, "", region, "Source term", varname, directdataname);
@@ -4989,13 +4987,11 @@ model_complex_plain_vector &
// cout << "is_lin : " << int(is_lin) << endl;
if (is_lin) {
- return add_linear_term
- (md, mim, expr, region, false, false,
- "Dirichlet condition with Nitsche's method");
+ return add_linear_term(md, mim, expr, region, false, false,
+ "Dirichlet condition with Nitsche's method");
} else {
- return add_nonlinear_term
- (md, mim, expr, region, false, false,
- "Dirichlet condition with Nitsche's method");
+ return add_nonlinear_term(md, mim, expr, region, false, false,
+ "Dirichlet condition with Nitsche's method");
}
}
@@ -5024,13 +5020,11 @@ model_complex_plain_vector &
+derivative_Neumann+")";
}
if (is_lin) {
- return add_linear_term
- (md, mim, expr, region, false, false,
- "Dirichlet condition with Nitsche's method");
+ return add_linear_term(md, mim, expr, region, false, false,
+ "Dirichlet condition with Nitsche's method");
} else {
- return add_nonlinear_term
- (md, mim, expr, region, false, false,
- "Dirichlet condition with Nitsche's method");
+ return add_nonlinear_term(md, mim, expr, region, false, false,
+ "Dirichlet condition with Nitsche's method");
}
}
@@ -5059,13 +5053,11 @@ model_complex_plain_vector &
+derivative_Neumann+"))";
}
if (is_lin) {
- return add_linear_term
- (md, mim, expr, region, false, false,
- "Dirichlet condition with Nitsche's method");
+ return add_linear_term(md, mim, expr, region, false, false,
+ "Dirichlet condition with Nitsche's method");
} else {
- return add_nonlinear_term
- (md, mim, expr, region, false, false,
- "Dirichlet condition with Nitsche's method");
+ return add_nonlinear_term(md, mim, expr, region, false, false,
+ "Dirichlet condition with Nitsche's method");
}
}
@@ -5463,11 +5455,11 @@ model_complex_plain_vector &
std::string expr = "Grad_"+varname+".Grad_"+test_varname
+" + sqr("+dataexpr+")*"+varname+"*"+test_varname;
- size_type ib = add_linear_term
- (md, mim, expr, region, true, true, "Helmholtz", true);
+ size_type ib = add_linear_term(md, mim, expr, region, true, true,
+ "Helmholtz", true);
if (ib == size_type(-1))
- ib = add_nonlinear_term
- (md, mim, expr, region, false, false, "Helmholtz (nonlinear)");
+ ib = add_nonlinear_term(md, mim, expr, region, false, false,
+ "Helmholtz (nonlinear)");
return ib;
}
}
@@ -5582,11 +5574,11 @@ model_complex_plain_vector &
std::string test_varname
= "Test_" + sup_previous_and_dot_to_varname(varname);
std::string expr = "(("+dataexpr+")*"+varname+")."+test_varname;
- size_type ib = add_linear_term
- (md, mim, expr, region, true, true, "Fourier-Robin", true);
+ size_type ib = add_linear_term(md, mim, expr, region, true, true,
+ "Fourier-Robin", true);
if (ib == size_type(-1))
- ib = add_nonlinear_term
- (md, mim, expr, region, false, false, "Fourier-Robin (nonlinear)");
+ ib = add_nonlinear_term(md, mim, expr, region, false, false,
+ "Fourier-Robin (nonlinear)");
return ib;
}
}
@@ -6088,8 +6080,8 @@ model_complex_plain_vector &
bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
if (is_lin) {
- return add_linear_term
- (md, mim, expr, region, false, false, "Linearized isotropic
elasticity");
+ return add_linear_term(md, mim, expr, region, false, false,
+ "Linearized isotropic elasticity");
} else {
return add_nonlinear_term
(md, mim, expr, region, false, false,
@@ -6122,8 +6114,8 @@ model_complex_plain_vector &
bool is_lin = workspace.used_variables(vl, vl_test1, vl_test2, dl, 2);
if (is_lin) {
- return add_linear_term
- (md, mim, expr, region, false, false, "Linearized isotropic
elasticity");
+ return add_linear_term(md, mim, expr, region, false, false,
+ "Linearized isotropic elasticity");
} else {
return add_nonlinear_term
(md, mim, expr, region, false, false,
@@ -6330,8 +6322,8 @@ model_complex_plain_vector &
else
expr = "-"+multname+"*Div_"+test_varname + "-"+test_multname
+"*Div_"+varname;
- size_type ib = add_linear_term
- (md, mim, expr, region, true, true, "Linear incompressibility", true);
+ size_type ib = add_linear_term(md, mim, expr, region, true, true,
+ "Linear incompressibility", true);
if (ib == size_type(-1))
ib = add_nonlinear_term
(md, mim, expr, region, false, false,
@@ -6471,11 +6463,11 @@ model_complex_plain_vector &
expr ="(("+dataexpr_rho+")*"+varname+")."+test_varname;
else
expr = varname+"."+test_varname;
- size_type ib = add_linear_term
- (md, mim, expr, region, true, true, "Mass matrix", true);
+ size_type ib = add_linear_term(md, mim, expr, region, true, true,
+ "Mass matrix", true);
if (ib == size_type(-1))
- ib = add_nonlinear_term
- (md, mim, expr, region, false, false, "Mass matrix (nonlinear)");
+ ib = add_nonlinear_term(md, mim, expr, region, false, false,
+ "Mass matrix (nonlinear)");
return ib;
}
}