branch: master commit a14b31c7f24724b7de18ba967f8bb129e18bec76 Author: Konstantinos Poulios <logar...@gmail.com> Date: Fri Nov 9 18:26:52 2018 +0100
Code clean up --- interface/src/gf_asm.cc | 2 +- src/bgeot_geometric_trans.cc | 20 ++++++------ src/getfem/getfem_projected_fem.h | 8 ++--- src/getfem_contact_and_friction_common.cc | 30 ++++++++--------- src/getfem_generic_assembly_semantic.cc | 19 +++++------ src/getfem_projected_fem.cc | 2 +- src/getfem_regular_meshes.cc | 53 ++++++++++++++----------------- 7 files changed, 64 insertions(+), 70 deletions(-) diff --git a/interface/src/gf_asm.cc b/interface/src/gf_asm.cc index d746be6..42416af 100644 --- a/interface/src/gf_asm.cc +++ b/interface/src/gf_asm.cc @@ -1254,7 +1254,7 @@ void gf_asm(getfemint::mexargs_in& m_in, getfemint::mexargs_out& m_out) { ); - /*@FUNC ('expression analysis', @str expression [, {@tm mesh | @tmim mim}] [, der_order] [, @tmodel model] [, @str varname, @int is_variable[, {@tmf mf | @tmimd mimd}], ...]) + /*@FUNC ('expression analysis', @str expression [, {@tmesh mesh | @tmim mim}] [, der_order] [, @tmodel model] [, @str varname, @int is_variable[, {@tmf mf | @tmimd mimd}], ...]) Analyse a high-level generic assembly expression and print information about the provided expression.@*/ sub_command diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc index 43d1c7d..7d29aa1 100644 --- a/src/bgeot_geometric_trans.cc +++ b/src/bgeot_geometric_trans.cc @@ -242,7 +242,7 @@ namespace bgeot { scalar_type a0 = A[4]*A[8] - A[5]*A[7], a1 = A[5]*A[6] - A[3]*A[8]; scalar_type a2 = A[3]*A[7] - A[4]*A[6]; scalar_type det = A[0] * a0 + A[1] * a1 + A[2] * a2; - GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix"); + GMM_ASSERT1(det != scalar_type(0), "Non invertible matrix"); scalar_type a3 = (A[2]*A[7] - A[1]*A[8]), a6 = (A[1]*A[5] - A[2]*A[4]); scalar_type a4 = (A[0]*A[8] - A[2]*A[6]), a7 = (A[2]*A[3] - A[0]*A[5]); scalar_type a5 = (A[1]*A[6] - A[0]*A[7]), a8 = (A[0]*A[4] - A[1]*A[3]); @@ -343,9 +343,9 @@ namespace bgeot { J__ = it[0] * a0 + it[1] * a1 + it[2] * a2; } break; default: - B_factors.base_resize(P, P); // store factorization for B computation - gmm::copy(gmm::transposed(KK), B_factors); - ipvt.resize(P); + B_factors.base_resize(P, P); // store factorization for B computation + gmm::copy(gmm::transposed(KK), B_factors); + ipvt.resize(P); bgeot::lu_factor(&(*(B_factors.begin())), ipvt, P); J__ = bgeot::lu_det(&(*(B_factors.begin())), ipvt, P); break; @@ -387,8 +387,10 @@ namespace bgeot { case 2: { auto it = &(*(KK.begin())); auto itB = &(*(B_.begin())); - *itB++ = it[3] / J__; *itB++ = -it[2] / J__; - *itB++ = -it[1] / J__; *itB = (*it) / J__; + *itB++ = it[3] / J__; + *itB++ = -it[2] / J__; + *itB++ = -it[1] / J__; + *itB = (*it) / J__; } break; case 3: { @@ -1073,11 +1075,9 @@ namespace bgeot { /* norm of returned vector is the ratio between the face surface on the reference element and the face surface on the real element. IT IS NOT UNITARY - - pt is the position of the evaluation point on the reference element */ - base_small_vector compute_normal(const geotrans_interpolation_context& c, - size_type face) { + base_small_vector + compute_normal(const geotrans_interpolation_context& c, size_type face) { GMM_ASSERT1(c.G().ncols() == c.pgt()->nb_points(), "dimensions mismatch"); base_small_vector un(c.N()); gmm::mult(c.B(), c.pgt()->normals()[face], un); diff --git a/src/getfem/getfem_projected_fem.h b/src/getfem/getfem_projected_fem.h index dc06452..fb8b3d6 100644 --- a/src/getfem/getfem_projected_fem.h +++ b/src/getfem/getfem_projected_fem.h @@ -95,8 +95,8 @@ namespace getfem { mutable bgeot::kdtree tree; // Tree containing the nodes of the // projected mf_source dofs mutable std::vector<size_type> ind_dof; /* all functions using this work - array should keep it full of - size_type(-1) */ + array should keep it filled + with size_type(-1) */ mutable bgeot::geotrans_inv_convex gic; mutable base_tensor taux; mutable fem_interpolation_context fictx; @@ -106,9 +106,9 @@ namespace getfem { mutable bgeot::multi_index mi2, mi3; mutable base_node ptref; - void build_kdtree(void) const; + void build_kdtree() const; - bool find_a_projected_point(base_node pt, base_node &ptr_proj, + bool find_a_projected_point(const base_node &pt, base_node &ptr_proj, size_type &cv_proj, short_type &fc_proj) const; virtual void update_from_context(void) const; diff --git a/src/getfem_contact_and_friction_common.cc b/src/getfem_contact_and_friction_common.cc index fada408..e4d05a6 100644 --- a/src/getfem_contact_and_friction_common.cc +++ b/src/getfem_contact_and_friction_common.cc @@ -1933,7 +1933,7 @@ namespace getfem { gmm::add(gmm::identity_matrix(), F_x); gmm::copy(F_x, F_x_inv); bgeot::lu_inverse(&(*(F_x_inv.begin())), N); - + base_tensor base_ux; base_matrix vbase_ux; @@ -1941,7 +1941,7 @@ namespace getfem { size_type qdim_ux = pfu_x->target_dim(); size_type ndof_ux = pfu_x->nb_dof(cv_x) * N / qdim_ux; vectorize_base_tensor(base_ux, vbase_ux, ndof_ux, qdim_ux, N); - + base_tensor base_uy; base_matrix vbase_uy; size_type ndof_uy = 0; @@ -1951,7 +1951,7 @@ namespace getfem { ndof_uy = pfu_y->nb_dof(cv_y) * N / qdim_uy; vectorize_base_tensor(base_uy, vbase_uy, ndof_uy, qdim_uy, N); } - + base_tensor grad_base_ux, vgrad_base_ux; ctx_x.grad_base_value(grad_base_ux); vectorize_grad_base_tensor(grad_base_ux, vgrad_base_ux, ndof_ux, @@ -1965,7 +1965,7 @@ namespace getfem { gmm::mult(F_y_inv, I_nxny, M1); base_matrix der_x(ndof_ux, N); gmm::mult(vbase_ux, gmm::transposed(M1), der_x); - + // -F_y^{-1}*I_nxny*Test_u(Y) base_matrix der_y(ndof_uy, N); if (ret_type == 1) { @@ -2743,22 +2743,22 @@ namespace getfem { scalar_type norm(0); if (tau > scalar_type(0)) { - gmm::add(lambda, gmm::scaled(Vs, -r), F); - scalar_type mu = gmm::vect_sp(F, n)/nn; - gmm::add(gmm::scaled(n, -mu/nn), F); + gmm::add(lambda, gmm::scaled(Vs, -r), F); // F <-- lambda -r*Vs + scalar_type mu = gmm::vect_sp(F, n)/nn; // mu <-- (lambda -r*Vs).n/|n| + gmm::add(gmm::scaled(n, -mu/nn), F); // F <-- (lambda -r*Vs)*(I-n x n / |n|²) norm = gmm::vect_norm2(F); - gmm::copy(gmm::identity_matrix(), dn); - gmm::scale(dn, -mu/nn); - gmm::rank_one_update(dn, gmm::scaled(n, mu/(nn*nn*nn)), n); - gmm::rank_one_update(dn, gmm::scaled(n, scalar_type(-1)/(nn*nn)), F); - gmm::copy(gmm::identity_matrix(), dVs); - gmm::rank_one_update(dVs, n, gmm::scaled(n, scalar_type(-1)/(nn*nn))); + gmm::copy(gmm::identity_matrix(), dn); // dn <-- I + gmm::scale(dn, -mu/nn); // dn <-- -(lambda -r*Vs).n/|n|² I + gmm::rank_one_update(dn, gmm::scaled(n, mu/(nn*nn*nn)), n); // dn <-- -(lambda -r*Vs).n/|n|² (I - n x n/|n|²) + gmm::rank_one_update(dn, gmm::scaled(n, scalar_type(-1)/(nn*nn)), F); // dn <-- -(lambda -r*Vs).n/|n|² (I - n x n/|n|²) + n x ((lambda -r*Vs)*(I-n x n / |n|²)) /|n|² + gmm::copy(gmm::identity_matrix(), dVs); // dVs <-- I + gmm::rank_one_update(dVs, n, gmm::scaled(n, scalar_type(-1)/(nn*nn))); // dVs <-- I - n x n/|n|² - if (norm > tau) { + if (norm > tau) { // slip gmm::rank_one_update(dVs, F, gmm::scaled(F, scalar_type(-1)/(norm*norm))); gmm::scale(dVs, tau / norm); - gmm::copy(gmm::scaled(F, scalar_type(1)/norm), dg); + gmm::copy(gmm::scaled(F, scalar_type(1)/norm), dg); // dg <-- Normalized((lambda -r*Vs)*(I-n x n / |n|²)) gmm::rank_one_update(dn, gmm::scaled(F, mu/(norm*norm*nn)), F); gmm::scale(dn, tau / norm); gmm::scale(F, tau / norm); diff --git a/src/getfem_generic_assembly_semantic.cc b/src/getfem_generic_assembly_semantic.cc index 524e05e..956f801 100644 --- a/src/getfem_generic_assembly_semantic.cc +++ b/src/getfem_generic_assembly_semantic.cc @@ -375,13 +375,13 @@ namespace getfem { } # define ga_valid_operand(pnode) \ - { \ - if (pnode && (pnode->node_type == GA_NODE_PREDEF_FUNC || \ - pnode->node_type == GA_NODE_SPEC_FUNC || \ - pnode->node_type == GA_NODE_NAME || \ - pnode->node_type == GA_NODE_OPERATOR || \ - pnode->node_type == GA_NODE_ALLINDICES)) \ - ga_throw_error(pnode->expr, pnode->pos, "Invalid term"); \ + { \ + if (pnode && (pnode->node_type == GA_NODE_PREDEF_FUNC || \ + pnode->node_type == GA_NODE_SPEC_FUNC || \ + pnode->node_type == GA_NODE_NAME || \ + pnode->node_type == GA_NODE_OPERATOR || \ + pnode->node_type == GA_NODE_ALLINDICES)) \ + ga_throw_error(pnode->expr, pnode->pos, "Invalid term"); \ } static void ga_node_analysis(ga_tree &tree, @@ -3175,7 +3175,7 @@ namespace getfem { pnode_trans = pnode->parent->children[1]; } - if (ivar) { + if (ivar) { // Derivative wrt the interpolated variable mi.resize(1); mi[0] = 2; for (size_type i = 0; i < pnode->tensor_order(); ++i) mi.push_back(pnode->tensor_proper_size(i)); @@ -3191,7 +3191,8 @@ namespace getfem { pnode->test_function_type = order; } - if (itrans) { + if (itrans) { // Derivative with respect to a variable that the + // interpolate transformation depends on const mesh_fem *mf = workspace.associated_mf(pnode_trans->name); size_type q = workspace.qdim(pnode_trans->name); size_type n = mf->linked_mesh().dim(); diff --git a/src/getfem_projected_fem.cc b/src/getfem_projected_fem.cc index 34adf8d..0050b0b 100644 --- a/src/getfem_projected_fem.cc +++ b/src/getfem_projected_fem.cc @@ -251,7 +251,7 @@ namespace getfem { tree.add_point_with_id(mf_source.point_of_basic_dof(dof), dof); } - bool projected_fem::find_a_projected_point(base_node pt, base_node &ptr_proj, + bool projected_fem::find_a_projected_point(const base_node &pt, base_node &ptr_proj, size_type &cv_proj, short_type &fc_proj) const { bgeot::index_node_pair ipt; diff --git a/src/getfem_regular_meshes.cc b/src/getfem_regular_meshes.cc index c319277..d856c18 100644 --- a/src/getfem_regular_meshes.cc +++ b/src/getfem_regular_meshes.cc @@ -294,17 +294,15 @@ namespace getfem std::string GT = PARAM.string_value("GT"); GMM_ASSERT1(!GT.empty(), "regular mesh : you have at least to " "specify the geometric transformation"); - bgeot::pgeometric_trans pgt = - bgeot::geometric_trans_descriptor(GT); + bgeot::pgeometric_trans pgt = bgeot::geometric_trans_descriptor(GT); size_type N = pgt->dim(); base_small_vector org(N); gmm::clear(org); - const std::vector<bgeot::md_param::param_value> &o - = PARAM.array_value("ORG"); + const auto &o = PARAM.array_value("ORG"); if (o.size() > 0) { - GMM_ASSERT1(o.size() == N, "ORG parameter should be an array of size " - << N); + GMM_ASSERT1(o.size() == N, + "ORG parameter should be an array of size " << N); for (size_type i = 0; i < N; ++i) { GMM_ASSERT1(o[i].type_of_param() == bgeot::md_param::REAL_VALUE, "ORG should be a real array."); @@ -316,14 +314,13 @@ namespace getfem std::vector<size_type> nsubdiv(N); gmm::fill(nsubdiv, 2); - const std::vector<bgeot::md_param::param_value> &ns - = PARAM.array_value("NSUBDIV"); + const auto &ns = PARAM.array_value("NSUBDIV"); if (ns.size() > 0) { GMM_ASSERT1(ns.size() == N, "NSUBDIV parameter should be an array of size " << N); for (size_type i = 0; i < N; ++i) { GMM_ASSERT1(ns[i].type_of_param() == bgeot::md_param::REAL_VALUE, - "NSUBDIV should be an integer array."); + "NSUBDIV should be an integer array"); nsubdiv[i] = size_type(ns[i].real()+0.5); } } @@ -331,14 +328,13 @@ namespace getfem base_small_vector sizes(N); gmm::fill(sizes, 1.0); - const std::vector<bgeot::md_param::param_value> &si - = PARAM.array_value("SIZES"); + const auto &si = PARAM.array_value("SIZES"); if (si.size() > 0) { GMM_ASSERT1(si.size() == N, "SIZES parameter should be an array of size " << N); for (size_type i = 0; i < N; ++i) { GMM_ASSERT1(si[i].type_of_param() == bgeot::md_param::REAL_VALUE, - "SIZES should be a real array."); + "SIZES should be a real array"); sizes[i] = si[i].real(); } } @@ -361,20 +357,18 @@ namespace getfem std::string GT = PARAM.string_value("GT"); GMM_ASSERT1(!GT.empty(), "regular ball mesh : you have at least to " "specify the geometric transformation"); - bgeot::pgeometric_trans pgt = - bgeot::geometric_trans_descriptor(GT); + bgeot::pgeometric_trans pgt = bgeot::geometric_trans_descriptor(GT); size_type N = pgt->dim(); base_small_vector org(N); - const std::vector<bgeot::md_param::param_value> &o - = PARAM.array_value("ORG"); + const auto &o = PARAM.array_value("ORG"); if (o.size() > 0) { - GMM_ASSERT1(o.size() == N, "ORG parameter should be an array of size " - << N); + GMM_ASSERT1(o.size() == N, + "ORG parameter should be an array of size " << N); for (size_type i = 0; i < N; ++i) { GMM_ASSERT1(o[i].type_of_param() == bgeot::md_param::REAL_VALUE, - "ORG should be a real array."); + "ORG should be a real array"); org[i] = o[i].real(); } } @@ -383,29 +377,26 @@ namespace getfem bool noised = (PARAM.int_value("NOISED") != 0); size_type nsubdiv0(2), nsubdiv1(2); - const std::vector<bgeot::md_param::param_value> &ns - = PARAM.array_value("NSUBDIV"); + const auto &ns = PARAM.array_value("NSUBDIV"); if (ns.size() > 0) { GMM_ASSERT1(ns.size() == 2, - "NSUBDIV parameter should be an array of size " << 2); + "NSUBDIV parameter should be an array of size 2"); for (size_type i = 0; i < 2; ++i) GMM_ASSERT1(ns[i].type_of_param() == bgeot::md_param::REAL_VALUE, - "NSUBDIV should be an integer array."); + "NSUBDIV should be an integer array"); nsubdiv0 = size_type(ns[0].real()+0.5); nsubdiv1 = size_type(ns[1].real()+0.5); } scalar_type radius(1), core_ratio(M_SQRT1_2); - const std::vector<bgeot::md_param::param_value> &si - = PARAM.array_value("SIZES"); + const auto &si = PARAM.array_value("SIZES"); if (si.size() > 0) { GMM_ASSERT1(si.size() == 1, - "SIZES parameter should be an array of size " << 1); + "SIZES parameter should be an array of size 1"); GMM_ASSERT1(si[0].type_of_param() == bgeot::md_param::REAL_VALUE, - "SIZES should be a real array."); + "SIZES should be a real array"); radius = si[0].real(); } - std::vector<size_type> nsubdiv(N); gmm::fill(nsubdiv, nsubdiv0); @@ -442,7 +433,8 @@ namespace getfem trsl[i] = core_ratio; mm[i].translation(trsl); for (dal::bv_visitor cv(mm[i].convex_index()); !cv.finished(); ++cv) - m.add_convex_by_points(mm[i].trans_of_convex(cv), mm[i].points_of_convex(cv).begin()); + m.add_convex_by_points(mm[i].trans_of_convex(cv), + mm[i].points_of_convex(cv).begin()); } std::vector<base_node> pts(m.points().card(), base_node(N)); @@ -503,7 +495,8 @@ namespace getfem m0.copy_from(m); m0.transformation(M); for (dal::bv_visitor cv(m0.convex_index()); !cv.finished(); ++cv) - m.add_convex_by_points(m0.trans_of_convex(cv), m0.points_of_convex(cv).begin()); + m.add_convex_by_points(m0.trans_of_convex(cv), + m0.points_of_convex(cv).begin()); } m.translation(org);