branch: master commit 75ef42f0da2be172db2f884b0afa9301ba6a43e9 Author: Konstantinos Poulios <logar...@gmail.com> AuthorDate: Fri Jan 3 23:45:46 2020 +0100
Basic infrastructure for defining internal variables - no condensation functionality yet --- src/getfem/getfem_generic_assembly.h | 14 +++-- src/getfem/getfem_models.h | 72 ++++++++++++++++--------- src/getfem_generic_assembly_compile_and_exec.cc | 12 +++-- src/getfem_generic_assembly_workspace.cc | 36 ++++++++++++- src/getfem_models.cc | 39 +++++++++----- 5 files changed, 125 insertions(+), 48 deletions(-) diff --git a/src/getfem/getfem_generic_assembly.h b/src/getfem/getfem_generic_assembly.h index 453e34e..9604225 100644 --- a/src/getfem/getfem_generic_assembly.h +++ b/src/getfem/getfem_generic_assembly.h @@ -265,7 +265,7 @@ namespace getfem { const model *md; const ga_workspace *parent_workspace; bool with_parent_variables; - size_type nb_prim_dof, nb_tmp_dof; + size_type nb_prim_dof, nb_intern_dof, first_intern_dof, nb_tmp_dof; void init(); @@ -280,6 +280,7 @@ namespace getfem { bgeot::multi_index qdims; // For data having a qdim different than the // qdim of the fem or im_data (dim per dof for // dof data) and for constant variables. + const bool is_internal; size_type qdim() const { size_type q = 1; @@ -289,9 +290,9 @@ namespace getfem { var_description(bool is_var, const mesh_fem *mf_, const im_data *imd_, gmm::sub_interval I_, const model_real_plain_vector *V_, - size_type Q) + size_type Q, bool is_intern_=false) : is_variable(is_var), is_fem_dofs(mf_ != 0), mf(mf_), imd(imd_), - I(I_), V(V_), qdims(1) + I(I_), V(V_), qdims(1), is_internal(is_intern_) { GMM_ASSERT1(Q > 0, "Bad dimension"); qdims[0] = Q; @@ -444,6 +445,9 @@ namespace getfem { void add_im_variable(const std::string &name, const im_data &imd, const gmm::sub_interval &I, const model_real_plain_vector &VV); + void add_internal_im_variable(const std::string &name, const im_data &imd, + const gmm::sub_interval &I, + const model_real_plain_vector &VV); void add_fixed_size_variable(const std::string &name, const gmm::sub_interval &I, const model_real_plain_vector &VV); @@ -463,6 +467,8 @@ namespace getfem { bool variable_exists(const std::string &name) const; + bool is_internal_variable(const std::string &name) const; + const std::string &variable_in_group(const std::string &group_name, const mesh &m) const; @@ -551,6 +557,8 @@ namespace getfem { bool include_empty_int_points() const; size_type nb_primary_dof() const { return nb_prim_dof; } + size_type nb_internal_dof() const { return nb_intern_dof; } + size_type first_internal_dof() const { return first_intern_dof; } size_type nb_temporary_dof() const { return nb_tmp_dof; } void add_temporary_interval_for_unreduced_variable(const std::string &name); diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h index 14decdb..5da43c6 100644 --- a/src/getfem/getfem_models.h +++ b/src/getfem/getfem_models.h @@ -151,6 +151,8 @@ namespace getfem { bool is_complex; // The variable is complex numbers bool is_affine_dependent; // The variable depends in an affine way // to another variable. + bool is_internal; // An internal variable defined on integration + // points, condensed out of the global system. bool is_fem_dofs; // The variable is the dofs of a fem size_type n_iter; // Number of versions of the variable stored. size_type n_temp_iter; // Number of additional temporary versions @@ -198,7 +200,7 @@ namespace getfem { const std::string &filter_var_ = std::string(""), mesh_im const *filter_mim_ = 0) : is_variable(is_var), is_disabled(false), is_complex(is_compl), - is_affine_dependent(false), + is_affine_dependent(false), is_internal(false), is_fem_dofs(mf_ != 0), n_iter(std::max(size_type(1), n_it)), n_temp_iter(0), default_iter(0), ptsc(0), @@ -525,6 +527,9 @@ namespace getfem { /** States if a name corresponds to a declared data. */ bool is_true_data(const std::string &name) const; + /** States if a variable is condensed out of the global system. */ + bool is_internal_variable(const std::string &name) const; + bool is_affine_dependent_variable(const std::string &name) const; const std::string &org_variable(const std::string &name) const; @@ -549,6 +554,13 @@ namespace getfem { the whole tangent system. */ bool is_coercive() const { return is_coercive_; } + /** Return true if the model has at least one internal variable. */ + bool has_internal_variables() const { + for (const auto &v : variables) + if (v.second.is_internal && !v.second.is_disabled) return true; + return false; + } + /** Return true if all the model terms do not affect the coercivity of the whole tangent system. */ bool is_symmetric() const { return is_symmetric_; } @@ -614,37 +626,41 @@ namespace getfem { model_complex_plain_vector & set_complex_constant_part(const std::string &name) const; + private: template<typename VECTOR, typename T> - void from_variables(VECTOR &V, T) const { + void from_variables(VECTOR &V, bool with_internal, T) const { for (const auto &v : variables) - if (v.second.is_variable && !(v.second.is_affine_dependent) - && !(v.second.is_disabled)) - gmm::copy(v.second.real_value[0], - gmm::sub_vector(V, v.second.I)); + if (v.second.is_variable && !v.second.is_affine_dependent + && !v.second.is_disabled + && (with_internal || !v.second.is_internal)) + gmm::copy(v.second.real_value[0], gmm::sub_vector(V, v.second.I)); } template<typename VECTOR, typename T> - void from_variables(VECTOR &V, std::complex<T>) const { + void from_variables(VECTOR &V, bool with_internal, std::complex<T>) const { for (const auto &v : variables) - if (v.second.is_variable && !(v.second.is_affine_dependent) - && !(v.second.is_disabled)) - gmm::copy(v.second.complex_value[0], - gmm::sub_vector(V, v.second.I)); + if (v.second.is_variable && !v.second.is_affine_dependent + && !v.second.is_disabled + && (with_internal || !v.second.is_internal)) + gmm::copy(v.second.complex_value[0], gmm::sub_vector(V, v.second.I)); } - template<typename VECTOR> void from_variables(VECTOR &V) const { + public: + template<typename VECTOR> + void from_variables(VECTOR &V, bool with_internal=false) const { typedef typename gmm::linalg_traits<VECTOR>::value_type T; context_check(); if (act_size_to_be_done) actualize_sizes(); - from_variables(V, T()); + from_variables(V, with_internal, T()); } + private: template<typename VECTOR, typename T> - void to_variables(const VECTOR &V, T) { + void to_variables(const VECTOR &V, bool with_internal, T) { for (auto &&v : variables) - if (v.second.is_variable && !(v.second.is_affine_dependent) - && !(v.second.is_disabled)) { - gmm::copy(gmm::sub_vector(V, v.second.I), - v.second.real_value[0]); + if (v.second.is_variable && !v.second.is_affine_dependent + && !v.second.is_disabled + && (with_internal || !v.second.is_internal)) { + gmm::copy(gmm::sub_vector(V, v.second.I), v.second.real_value[0]); v.second.v_num_data[0] = act_counter(); } update_affine_dependent_variables(); @@ -652,22 +668,24 @@ namespace getfem { } template<typename VECTOR, typename T> - void to_variables(const VECTOR &V, std::complex<T>) { + void to_variables(const VECTOR &V, bool with_internal, std::complex<T>) { for (auto &&v : variables) - if (v.second.is_variable && !(v.second.is_affine_dependent) - && !(v.second.is_disabled)) { - gmm::copy(gmm::sub_vector(V, v.second.I), - v.second.complex_value[0]); + if (v.second.is_variable && !v.second.is_affine_dependent + && !v.second.is_disabled + && (with_internal || !v.second.is_internal)) { + gmm::copy(gmm::sub_vector(V, v.second.I), v.second.complex_value[0]); v.second.v_num_data[0] = act_counter(); } update_affine_dependent_variables(); this->post_to_variables_step(); } - template<typename VECTOR> void to_variables(const VECTOR &V) { + public: + template<typename VECTOR> + void to_variables(const VECTOR &V, bool with_internal=false) { typedef typename gmm::linalg_traits<VECTOR>::value_type T; context_check(); if (act_size_to_be_done) actualize_sizes(); - to_variables(V, T()); + to_variables(V, with_internal, T()); } /** Add a fixed size variable to the model assumed to be a vector. @@ -753,7 +771,9 @@ namespace getfem { /** Add variable defined at integration points. */ void add_im_variable(const std::string &name, const im_data &im_data, size_type niter = 1); - + /** Add internal variable, defined at integration points and condensated. */ + void add_internal_im_variable(const std::string &name, + const im_data &im_data); /** Add data defined at integration points. */ void add_im_data(const std::string &name, const im_data &im_data, size_type niter = 1); diff --git a/src/getfem_generic_assembly_compile_and_exec.cc b/src/getfem_generic_assembly_compile_and_exec.cc index 9dda8b4..3bb6e25 100644 --- a/src/getfem_generic_assembly_compile_and_exec.cc +++ b/src/getfem_generic_assembly_compile_and_exec.cc @@ -7251,10 +7251,11 @@ namespace getfem { GMM_ASSERT1(root->interpolate_name_test1.size() == 0, "Interpolate transformation on integration " "point variable"); - pgai = std::make_shared<ga_instruction_vector_assembly_imd> - (root->tensor(), Vr, gis.ctx, - workspace.interval_of_variable(root->name_test1), - *imd, gis.coeff, gis.ipt); + if (!workspace.is_internal_variable(root->name_test1)) + pgai = std::make_shared<ga_instruction_vector_assembly_imd> + (root->tensor(), Vr, gis.ctx, + workspace.interval_of_variable(root->name_test1), + *imd, gis.coeff, gis.ipt); } else { pgai = std::make_shared<ga_instruction_vector_assembly> (root->tensor(), Vr, @@ -7341,7 +7342,8 @@ namespace getfem { <ga_instruction_matrix_assembly_standard_vector> (root->tensor(), Krr, ctx1, ctx2, I1, I2, mf1, mf2, alpha1, alpha2, gis.coeff, gis.nbpt, gis.ipt); - } else { + } else if (!workspace.is_internal_variable(root->name_test1) && + !workspace.is_internal_variable(root->name_test2)) { auto &Kxu = (mf1 && mf1->is_reduced()) ? Kuu : Kru; auto &Kxr = (mf1 && mf1->is_reduced()) ? Kur : Krr; auto &Kux = (mf2 && mf2->is_reduced()) ? Kuu : Kur; diff --git a/src/getfem_generic_assembly_workspace.cc b/src/getfem_generic_assembly_workspace.cc index 6cd17db..c3c2ae1 100644 --- a/src/getfem_generic_assembly_workspace.cc +++ b/src/getfem_generic_assembly_workspace.cc @@ -72,6 +72,8 @@ namespace getfem { void ga_workspace::add_fem_variable (const std::string &name, const mesh_fem &mf, const gmm::sub_interval &I, const model_real_plain_vector &VV) { + GMM_ASSERT1(nb_intern_dof == 0 || I.last() < first_intern_dof, + "The provided interval overlaps with internal dofs"); nb_prim_dof = std::max(nb_prim_dof, I.last()); variables.emplace(name, var_description(true, &mf, 0, I, &VV, 1)); } @@ -79,13 +81,29 @@ namespace getfem { void ga_workspace::add_im_variable (const std::string &name, const im_data &imd, const gmm::sub_interval &I, const model_real_plain_vector &VV) { + GMM_ASSERT1(nb_intern_dof == 0 || I.last() <= first_intern_dof, + "The provided interval overlaps with internal dofs"); nb_prim_dof = std::max(nb_prim_dof, I.last()); variables.emplace(name, var_description(true, 0, &imd, I, &VV, 1)); } + void ga_workspace::add_internal_im_variable + (const std::string &name, const im_data &imd, + const gmm::sub_interval &I, const model_real_plain_vector &VV) { + GMM_ASSERT1(I.first() >= nb_prim_dof, + "The provided interval overlaps with primary dofs"); + nb_intern_dof += first_intern_dof - std::min(first_intern_dof, I.first()); + first_intern_dof = std::min(first_intern_dof, I.first()); + nb_intern_dof += first_intern_dof + nb_intern_dof + - std::min(first_intern_dof + nb_intern_dof, I.last()); + variables.emplace(name, var_description(true, 0, &imd, I, &VV, 1, true)); + } + void ga_workspace::add_fixed_size_variable (const std::string &name, const gmm::sub_interval &I, const model_real_plain_vector &VV) { + GMM_ASSERT1(nb_intern_dof == 0 || I.last() <= first_intern_dof, + "The provided interval overlaps with internal dofs"); nb_prim_dof = std::max(nb_prim_dof, I.last()); variables.emplace(name, var_description(true, 0, 0, I, &VV, dim_type(gmm::vect_size(VV)))); @@ -116,6 +134,18 @@ namespace getfem { gmm::vect_size(VV)/(imd.nb_filtered_index() * imd.nb_tensor_elem()))); } + bool ga_workspace::is_internal_variable(const std::string &name) const { + + if ((md && md->variable_exists(name) && md->is_internal_variable(name)) || + (parent_workspace && parent_workspace->variable_exists(name) + && parent_workspace->is_internal_variable(name))) + return true; + else { + VAR_SET::const_iterator it = variables.find(name); + return it == variables.end() ? false : it->second.is_internal; + } + } + bool ga_workspace::variable_exists(const std::string &name) const { return (md && md->variable_exists(name)) || (parent_workspace && parent_workspace->variable_exists(name)) || @@ -942,6 +972,7 @@ namespace getfem { { init(); nb_prim_dof = with_parent_variables ? md->nb_dof() : 0; + nb_intern_dof = 0; if (var_inherit == inherit::ALL) { // enable model's disabled variables model::varnamelist vlmd; md->variable_list(vlmd); @@ -965,6 +996,7 @@ namespace getfem { } } } + first_intern_dof = nb_prim_dof; // dofs are contiguous in getfem::model } ga_workspace::ga_workspace(const ga_workspace &gaw, const inherit var_inherit) @@ -975,10 +1007,12 @@ namespace getfem { { init(); nb_prim_dof = with_parent_variables ? gaw.nb_primary_dof() : 0; + nb_intern_dof = with_parent_variables ? gaw.nb_internal_dof() : 0; + first_intern_dof = with_parent_variables ? gaw.first_internal_dof() : 0; } ga_workspace::ga_workspace() : md(0), parent_workspace(0), with_parent_variables(false), - nb_prim_dof(0), nb_tmp_dof(0) + nb_prim_dof(0), nb_intern_dof(0), first_intern_dof(0), nb_tmp_dof(0) { init(); } ga_workspace::~ga_workspace() { clear_expressions(); } diff --git a/src/getfem_models.cc b/src/getfem_models.cc index 725094b..d0330eb 100644 --- a/src/getfem_models.cc +++ b/src/getfem_models.cc @@ -244,6 +244,12 @@ namespace getfem { return is_old(name) || !(variable_description(name).is_variable); } + bool model::is_internal_variable(const std::string &name) const { + if (is_old(name)) return false; + const auto &var_descr = variable_description(name); + return var_descr.is_internal && !var_descr.is_disabled; + } + bool model::is_affine_dependent_variable(const std::string &name) const { return !(is_old(name)) && variable_description(name).is_affine_dependent; } @@ -293,17 +299,24 @@ namespace getfem { } void model::resize_global_system() const { - size_type tot_size = 0; - for (auto &&v : variables) { - if (v.second.is_variable && v.second.is_disabled) - v.second.I = gmm::sub_interval(0,0); - if (v.second.is_variable && !(v.second.is_affine_dependent) - && !(v.second.is_disabled)) { - v.second.I = gmm::sub_interval(tot_size, v.second.size()); - tot_size += v.second.size(); + size_type full_size = 0; + for (auto &&v : variables) + if (v.second.is_variable) { + if (v.second.is_disabled) + v.second.I = gmm::sub_interval(0,0); + else if (!v.second.is_affine_dependent && !v.second.is_internal) { + v.second.I = gmm::sub_interval(full_size, v.second.size()); + full_size += v.second.size(); + } + } + size_type primary_size = full_size; + + for (auto &&v : variables) + if (v.second.is_internal && !v.second.is_disabled) { // is_internal_variable() + v.second.I = gmm::sub_interval(full_size, v.second.size()); + full_size += v.second.size(); } - } for (auto &&v : variables) if (v.second.is_affine_dependent) { @@ -312,12 +325,12 @@ namespace getfem { } if (complex_version) { - gmm::resize(cTM, tot_size, tot_size); - gmm::resize(crhs, tot_size); + gmm::resize(cTM, primary_size, primary_size); + gmm::resize(crhs, primary_size); } else { - gmm::resize(rTM, tot_size, tot_size); - gmm::resize(rrhs, tot_size); + gmm::resize(rTM, primary_size, primary_size); + gmm::resize(rrhs, primary_size); } for (dal::bv_visitor ib(valid_bricks); !ib.finished(); ++ib)