[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sun, 29 Apr 2018 13:11:23 -0400 (EDT) |
branch: devel-yves-generic-assembly-modifs
commit ee8ed293a7262b3455f324fdc22bbadcb1bfe8c8
Author: Yves Renard <address@hidden>
Date: Sat Apr 28 15:51:32 2018 +0200
split of getfem_generic_assembly.cc done
---
src/Makefile.am | 4 +-
src/getfem/getfem_generic_assembly.h | 1 +
.../getfem_generic_assembly_compile_and_exec.h | 24 +-
...tfem_generic_assembly_functions_and_operators.h | 43 +-
src/getfem/getfem_generic_assembly_semantic.h | 2 +-
...=> getfem_generic_assembly_compile_and_exec.cc} | 1955 +-------------------
...fem_generic_assembly_functions_and_operators.cc | 115 +-
src/getfem_generic_assembly_interpolation.cc | 914 +++++++++
src/getfem_generic_assembly_semantic.cc | 17 +-
src/getfem_generic_assembly_tree.cc | 10 +-
src/getfem_generic_assembly_workspace.cc | 1005 ++++++++++
11 files changed, 2090 insertions(+), 2000 deletions(-)
diff --git a/src/Makefile.am b/src/Makefile.am
index 30ed222..4d14864 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -223,10 +223,12 @@ SRC =
\
getfem_error_estimate.cc \
getfem_export.cc \
getfem_assembling_tensors.cc \
- getfem_generic_assembly.cc \
getfem_generic_assembly_tree.cc \
getfem_generic_assembly_functions_and_operators.cc \
getfem_generic_assembly_semantic.cc \
+ getfem_generic_assembly_workspace.cc \
+ getfem_generic_assembly_compile_and_exec.cc \
+ getfem_generic_assembly_interpolation.cc \
getfem_mesher.cc \
getfem_fourth_order.cc \
getfem_nonlinear_elasticity.cc \
diff --git a/src/getfem/getfem_generic_assembly.h
b/src/getfem/getfem_generic_assembly.h
index 1d8d1b9..f7d229a 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -160,6 +160,7 @@ namespace getfem {
void add_method(const std::string &name,
const std::shared_ptr<ga_nonlinear_operator> &pt)
{ tab[name] = pt; }
+ ga_predef_operator_tab();
};
//=========================================================================
diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h
b/src/getfem/getfem_generic_assembly_compile_and_exec.h
index 66f8757..0f7766d 100644
--- a/src/getfem/getfem_generic_assembly_compile_and_exec.h
+++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h
@@ -82,20 +82,7 @@ namespace getfem {
};
bool operator <(const gauss_pt_corresp &gpc1,
- const gauss_pt_corresp &gpc2) {
- if (gpc1.pai != gpc2.pai)
- return (gpc1.pai < gpc2.pai );
- if (gpc1.nodes.size() != gpc2.nodes.size())
- return (gpc1.nodes.size() < gpc2.nodes.size());
- for (size_type i = 0; i < gpc1.nodes.size(); ++i)
- if (gpc1.nodes[i] != gpc2.nodes[i])
- return (gpc1.nodes[i] < gpc2.nodes[i]);
- if (gpc1.pgt1 != gpc2.pgt1)
- return (gpc1.pgt1 < gpc2.pgt1);
- if (gpc1.pgt2 != gpc2.pgt2)
- return (gpc1.pgt2 < gpc2.pgt2);
- return false;
- }
+ const gauss_pt_corresp &gpc2);
struct ga_instruction_set {
@@ -210,6 +197,15 @@ namespace getfem {
size_type order);
void ga_compile_function(ga_workspace &workspace,
ga_instruction_set &gis, bool scalar);
+ void ga_compile_interpolation(ga_workspace &workspace,
+ ga_instruction_set &gis);
+ void ga_interpolation_exec(ga_instruction_set &gis,
+ ga_workspace &workspace,
+ ga_interpolation_context &gic);
+ void ga_interpolation_single_point_exec
+ (ga_instruction_set &gis, ga_workspace &workspace,
+ const fem_interpolation_context &ctx_x, const base_small_vector &Normal,
+ const mesh &interp_mesh);
} /* end of namespace */
diff --git a/src/getfem/getfem_generic_assembly_functions_and_operators.h
b/src/getfem/getfem_generic_assembly_functions_and_operators.h
index 9039a7c..ca2f06a 100644
--- a/src/getfem/getfem_generic_assembly_functions_and_operators.h
+++ b/src/getfem/getfem_generic_assembly_functions_and_operators.h
@@ -78,32 +78,33 @@ namespace getfem {
bool is_affine(const std::string &varname) const;
- size_type ftype() const {return ftype_;}
- size_type dtype() const {return dtype_;}
- size_type nbargs() const {return nbargs_;}
- const std::string &derivative1() const {return derivative1_;}
- const std::string &derivative2() const {return derivative2_;}
- const std::string &expr() const {return expr_;}
- pscalar_func_onearg f1() const {return f1_;}
- pscalar_func_twoargs f2() const {return f2_;}
-
- ga_predef_function() : expr_(""), derivative1_(""),
- derivative2_(""), gis(nullptr) {}
+ size_type ftype() const { return ftype_;}
+ size_type dtype() const { return dtype_;}
+ size_type nbargs() const { return nbargs_;}
+ const std::string &derivative1() const { return derivative1_;}
+ const std::string &derivative2() const { return derivative2_;}
+ const std::string &expr() const { return expr_;}
+ pscalar_func_onearg f1() const { return f1_;}
+ pscalar_func_twoargs f2() const { return f2_;}
+
+ ga_predef_function();
ga_predef_function(pscalar_func_onearg f, size_type dtype__ = 0,
- const std::string &der = "")
- : ftype_(0), dtype_(dtype__), nbargs_(1), f1_(f), expr_(""),
- derivative1_(der), derivative2_("") {}
+ const std::string &der = "");
ga_predef_function(pscalar_func_twoargs f, size_type dtype__ = 0,
const std::string &der1 = "",
- const std::string &der2 = "")
- : ftype_(0), dtype_(dtype__), nbargs_(2), f2_(f),
- expr_(""), derivative1_(der1), derivative2_(der2), gis(nullptr) {}
- ga_predef_function(const std::string &expr__)
- : ftype_(1), dtype_(3), nbargs_(1), expr_(expr__),
- derivative1_(""), derivative2_(""), t(1, 0.), u(1, 0.), gis(nullptr) {}
+ const std::string &der2 = "");
+ ga_predef_function(const std::string &expr__);
};
- typedef std::map<std::string, ga_predef_function> ga_predef_function_tab;
+ struct ga_predef_function_tab
+ : public std::map<std::string, ga_predef_function> {
+
+ ga_predef_function_tab();
+ };
+
+ struct ga_spec_function_tab : public std::set<std::string> {
+ ga_spec_function_tab();
+ };
} /* end of namespace */
diff --git a/src/getfem/getfem_generic_assembly_semantic.h
b/src/getfem/getfem_generic_assembly_semantic.h
index 0c3939a..78abcfe 100644
--- a/src/getfem/getfem_generic_assembly_semantic.h
+++ b/src/getfem/getfem_generic_assembly_semantic.h
@@ -60,7 +60,7 @@ namespace getfem {
size_type meshdim,
size_type ref_elt_dim,
bool eval_fixed_size,
- bool ignore_X, int option);
+ bool ignore_X, int option = 0);
/* Extract the variables used in a sub-tree, ignoring or not the data.
Variable groups are taken into account. Return if at least one variable
diff --git a/src/getfem_generic_assembly.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
similarity index 77%
rename from src/getfem_generic_assembly.cc
rename to src/getfem_generic_assembly_compile_and_exec.cc
index 41ba5de..507db9e 100644
--- a/src/getfem_generic_assembly.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -22,13 +22,27 @@
#include "getfem/getfem_generic_assembly_tree.h"
#include "getfem/getfem_generic_assembly_semantic.h"
#include "getfem/getfem_generic_assembly_compile_and_exec.h"
+#include "getfem/getfem_generic_assembly_functions_and_operators.h"
namespace getfem {
- extern bool predef_operators_nonlinear_elasticity_initialized;
- extern bool predef_operators_plasticity_initialized;
- extern bool predef_operators_contact_initialized;
+ bool operator <(const gauss_pt_corresp &gpc1,
+ const gauss_pt_corresp &gpc2) {
+ if (gpc1.pai != gpc2.pai)
+ return (gpc1.pai < gpc2.pai );
+ if (gpc1.nodes.size() != gpc2.nodes.size())
+ return (gpc1.nodes.size() < gpc2.nodes.size());
+ for (size_type i = 0; i < gpc1.nodes.size(); ++i)
+ if (gpc1.nodes[i] != gpc2.nodes[i])
+ return (gpc1.nodes[i] < gpc2.nodes[i]);
+ if (gpc1.pgt1 != gpc2.pgt1)
+ return (gpc1.pgt1 < gpc2.pgt1);
+ if (gpc1.pgt2 != gpc2.pgt2)
+ return (gpc1.pgt2 < gpc2.pgt2);
+ return false;
+ }
+
//=========================================================================
// Instructions for compilation: basic optimized operations on tensors
//=========================================================================
@@ -4298,982 +4312,8 @@ namespace getfem {
- //=========================================================================
- // Structure dealing with user defined environment : constant, variables,
- // functions, operators.
- //=========================================================================
-
- void ga_workspace::init() {
- // allocate own storage for K an V to be used unless/until external
- // storage is provided with set_assembled_matrix/vector
- K = std::make_shared<model_real_sparse_matrix>(2,2);
- V = std::make_shared<base_vector>(2);
- // add default transformations
- add_interpolate_transformation
- ("neighbour_elt", interpolate_transformation_neighbour_instance());
- }
-
- // variables and variable groups
- 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) {
- variables[name] = var_description(true, true, &mf, I, &VV, 0, 1);
- }
-
- void ga_workspace::add_fixed_size_variable
- (const std::string &name,
- const gmm::sub_interval &I, const model_real_plain_vector &VV) {
- variables[name] = var_description(true, false, 0, I, &VV, 0,
- dim_type(gmm::vect_size(VV)));
- }
-
- void ga_workspace::add_fem_constant
- (const std::string &name, const mesh_fem &mf,
- const model_real_plain_vector &VV) {
- GMM_ASSERT1(mf.nb_dof(), "The provided mesh_fem of variable" << name
- << "has zero degrees of freedom.");
- size_type Q = gmm::vect_size(VV)/mf.nb_dof();
- if (Q == 0) Q = size_type(1);
- variables[name] = var_description(false, true, &mf,
- gmm::sub_interval(), &VV, 0, Q);
- }
-
- void ga_workspace::add_fixed_size_constant
- (const std::string &name, const model_real_plain_vector &VV) {
- variables[name] = var_description(false, false, 0,
- gmm::sub_interval(), &VV, 0,
- gmm::vect_size(VV));
- }
-
- void ga_workspace::add_im_data(const std::string &name, const im_data &imd,
- const model_real_plain_vector &VV) {
- variables[name] = var_description
- (false, false, 0, gmm::sub_interval(), &VV, &imd,
- gmm::vect_size(VV)/(imd.nb_filtered_index() * imd.nb_tensor_elem()));
- }
-
-
- bool ga_workspace::variable_exists(const std::string &name) const {
- return (md && md->variable_exists(name)) ||
- (parent_workspace && parent_workspace->variable_exists(name)) ||
- (variables.find(name) != variables.end());
- }
-
- bool ga_workspace::variable_group_exists(const std::string &name) const {
- return (variable_groups.find(name) != variable_groups.end()) ||
- (md && md->variable_group_exists(name)) ||
- (parent_workspace && parent_workspace->variable_group_exists(name));
- }
-
- const std::vector<std::string>&
- ga_workspace::variable_group(const std::string &group_name) const {
- std::map<std::string, std::vector<std::string> >::const_iterator
- it = variable_groups.find(group_name);
- if (it != variable_groups.end())
- return (variable_groups.find(group_name))->second;
- if (md && md->variable_group_exists(group_name))
- return md->variable_group(group_name);
- if (parent_workspace &&
- parent_workspace->variable_group_exists(group_name))
- return parent_workspace->variable_group(group_name);
- GMM_ASSERT1(false, "Undefined variable group " << group_name);
- }
-
- const std::string&
- ga_workspace::first_variable_of_group(const std::string &name) const {
- const std::vector<std::string> &t = variable_group(name);
- GMM_ASSERT1(t.size(), "Variable group " << name << " is empty");
- return t[0];
- }
-
- bool ga_workspace::is_constant(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) return !(it->second.is_variable);
- if (variable_group_exists(name))
- return is_constant(first_variable_of_group(name));
- if (md && md->variable_exists(name)) {
- if (enable_all_md_variables) return md->is_true_data(name);
- return md->is_data(name);
- }
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->is_constant(name);
- GMM_ASSERT1(false, "Undefined variable " << name);
- }
-
- bool ga_workspace::is_disabled_variable(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) return false;
- if (variable_group_exists(name))
- return is_disabled_variable(first_variable_of_group(name));
- if (md && md->variable_exists(name)) {
- if (enable_all_md_variables) return false;
- return md->is_disabled_variable(name);
- }
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->is_disabled_variable(name);
- GMM_ASSERT1(false, "Undefined variable " << name);
- }
-
- const scalar_type &
- ga_workspace::factor_of_variable(const std::string &name) const {
- static const scalar_type one(1);
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) return one;
- if (variable_group_exists(name))
- return one;
- if (md && md->variable_exists(name)) return md->factor_of_variable(name);
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->factor_of_variable(name);
- GMM_ASSERT1(false, "Undefined variable " << name);
- }
-
- const gmm::sub_interval &
- ga_workspace::interval_of_disabled_variable(const std::string &name) const {
- std::map<std::string, gmm::sub_interval>::const_iterator
- it1 = int_disabled_variables.find(name);
- if (it1 != int_disabled_variables.end()) return it1->second;
- if (md->is_affine_dependent_variable(name))
- return interval_of_disabled_variable(md->org_variable(name));
-
- size_type first = md->nb_dof();
- for (const std::pair<std::string, gmm::sub_interval> &p
- : int_disabled_variables)
- first = std::max(first, p.second.last());
-
- int_disabled_variables[name]
- = gmm::sub_interval(first, gmm::vect_size(value(name)));
- return int_disabled_variables[name];
- }
-
- const gmm::sub_interval &
- ga_workspace::interval_of_variable(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) return it->second.I;
- if (md && md->variable_exists(name)) {
- if (enable_all_md_variables && md->is_disabled_variable(name))
- return interval_of_disabled_variable(name);
- return md->interval_of_variable(name);
- }
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->interval_of_variable(name);
- GMM_ASSERT1(false, "Undefined variable " << name);
- }
-
- const mesh_fem *
- ga_workspace::associated_mf(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end())
- return it->second.is_fem_dofs ? it->second.mf : 0;
- if (md && md->variable_exists(name))
- return md->pmesh_fem_of_variable(name);
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->associated_mf(name);
- if (variable_group_exists(name))
- return associated_mf(first_variable_of_group(name));
- GMM_ASSERT1(false, "Undefined variable or group " << name);
- }
-
- const im_data *
- ga_workspace::associated_im_data(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) return it->second.imd;
- if (md && md->variable_exists(name))
- return md->pim_data_of_variable(name);
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->associated_im_data(name);
- if (variable_group_exists(name)) return 0;
- GMM_ASSERT1(false, "Undefined variable " << name);
- }
-
- size_type ga_workspace::qdim(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) {
- const mesh_fem *mf = it->second.is_fem_dofs ? it->second.mf : 0;
- const im_data *imd = it->second.imd;
- size_type n = it->second.qdim();
- if (mf) {
- return n * mf->get_qdim();
- } else if (imd) {
- return n * imd->tensor_size().total_size();
- }
- return n;
- }
- if (md && md->variable_exists(name))
- return md->qdim_of_variable(name);
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->qdim(name);
- if (variable_group_exists(name))
- return qdim(first_variable_of_group(name));
- GMM_ASSERT1(false, "Undefined variable or group " << name);
- }
-
- bgeot::multi_index
- ga_workspace::qdims(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end()) {
- const mesh_fem *mf = it->second.is_fem_dofs ? it->second.mf : 0;
- const im_data *imd = it->second.imd;
- size_type n = it->second.qdim();
- if (mf) {
- bgeot::multi_index mi = mf->get_qdims();
- if (n > 1 || it->second.qdims.size() > 1) {
- size_type i = 0;
- if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
- for (; i < it->second.qdims.size(); ++i)
- mi.push_back(it->second.qdims[i]);
- }
- return mi;
- } else if (imd) {
- bgeot::multi_index mi = imd->tensor_size();
- size_type q = n / imd->nb_filtered_index();
- GMM_ASSERT1(q % imd->nb_tensor_elem() == 0,
- "Invalid mesh im data vector");
- if (n > 1 || it->second.qdims.size() > 1) {
- size_type i = 0;
- if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
- for (; i < it->second.qdims.size(); ++i)
- mi.push_back(it->second.qdims[i]);
- }
- return mi;
- }
- return it->second.qdims;
- }
- if (md && md->variable_exists(name))
- return md->qdims_of_variable(name);
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->qdims(name);
- if (variable_group_exists(name))
- return qdims(first_variable_of_group(name));
- GMM_ASSERT1(false, "Undefined variable or group " << name);
- }
-
- const model_real_plain_vector &
- ga_workspace::value(const std::string &name) const {
- VAR_SET::const_iterator it = variables.find(name);
- if (it != variables.end())
- return *(it->second.V);
- if (md && md->variable_exists(name))
- return md->real_variable(name);
- if (parent_workspace && parent_workspace->variable_exists(name))
- return parent_workspace->value(name);
- if (variable_group_exists(name))
- return value(first_variable_of_group(name));
- GMM_ASSERT1(false, "Undefined variable or group " << name);
- }
-
- scalar_type ga_workspace::get_time_step() const {
- if (md) return md->get_time_step();
- if (parent_workspace) return parent_workspace->get_time_step();
- GMM_ASSERT1(false, "No time step defined here");
- }
- // Macros
- bool ga_workspace::macro_exists(const std::string &name) const {
- if (macros.find(name) != macros.end()) return true;
- if (md && md->macro_exists(name)) return true;
- if (parent_workspace &&
- parent_workspace->macro_exists(name)) return true;
- return false;
- }
-
- const std::string&
- ga_workspace::get_macro(const std::string &name) const {
- std::map<std::string, std::string>::const_iterator it=macros.find(name);
- if (it != macros.end()) return it->second;
- if (md && md->macro_exists(name)) return md->get_macro(name);
- if (parent_workspace &&
- parent_workspace->macro_exists(name))
- return parent_workspace->get_macro(name);
- GMM_ASSERT1(false, "Undefined macro");
- }
-
- ga_tree &
- ga_workspace::macro_tree(const std::string &name, size_type meshdim,
- size_type ref_elt_dim, bool ignore_X) const {
- GMM_ASSERT1(macro_exists(name), "Undefined macro");
- auto it = macro_trees.find(name);
- bool to_be_analyzed = false;
- m_tree *mt = 0;
-
- if (it == macro_trees.end()) {
- mt = &(macro_trees[name]);
- to_be_analyzed = true;
- } else {
- mt = &(it->second);
- GMM_ASSERT1(mt->ptree, "Recursive definition of macro " << name);
- if (mt->meshdim != meshdim || mt->ignore_X != ignore_X) {
- to_be_analyzed = true;
- delete mt->ptree; mt->ptree = 0;
- }
- }
- if (to_be_analyzed) {
- ga_tree tree;
- ga_read_string(get_macro(name), tree);
- ga_semantic_analysis(get_macro(name), tree, *this, meshdim, ref_elt_dim,
- false, ignore_X, 3);
- GMM_ASSERT1(tree.root, "Invalid macro");
- mt->ptree = new ga_tree(tree);
- mt->meshdim = meshdim;
- mt->ignore_X = ignore_X;
- }
- return *(mt->ptree);
- }
-
- void ga_workspace::add_interpolate_transformation
- (const std::string &name, pinterpolate_transformation ptrans) {
- if (transformations.find(name) != transformations.end())
- GMM_ASSERT1(name.compare("neighbour_elt"), "neighbour_elt is a "
- "reserved interpolate transformation name");
- transformations[name] = ptrans;
- }
-
- bool ga_workspace::interpolate_transformation_exists
- (const std::string &name) const {
- return (md && md->interpolate_transformation_exists(name)) ||
- (parent_workspace &&
- parent_workspace->interpolate_transformation_exists(name)) ||
- (transformations.find(name) != transformations.end());
- }
-
- pinterpolate_transformation
- ga_workspace::interpolate_transformation(const std::string &name) const {
- std::map<std::string, pinterpolate_transformation>::const_iterator
- it = transformations.find(name);
- if (it != transformations.end()) return it->second;
- if (md && md->interpolate_transformation_exists(name))
- return md->interpolate_transformation(name);
- if (parent_workspace &&
- parent_workspace->interpolate_transformation_exists(name))
- return parent_workspace->interpolate_transformation(name);
- GMM_ASSERT1(false, "Inexistent transformation " << name);
- }
-
- bool ga_workspace::elementary_transformation_exists
- (const std::string &name) const {
- return (md && md->elementary_transformation_exists(name)) ||
- (parent_workspace &&
- parent_workspace->elementary_transformation_exists(name)) ||
- (elem_transformations.find(name) != elem_transformations.end());
- }
-
- pelementary_transformation
- ga_workspace::elementary_transformation(const std::string &name) const {
- std::map<std::string, pelementary_transformation>::const_iterator
- it = elem_transformations.find(name);
- if (it != elem_transformations.end()) return it->second;
- if (md && md->elementary_transformation_exists(name))
- return md->elementary_transformation(name);
- if (parent_workspace &&
- parent_workspace->elementary_transformation_exists(name))
- return parent_workspace->elementary_transformation(name);
- GMM_ASSERT1(false, "Inexistent elementary transformation " << name);
- }
-
- const mesh_region &
- ga_workspace::register_region(const mesh &m, const mesh_region ®ion) {
- if (&m == &dummy_mesh())
- return dummy_mesh_region();
-
- std::list<mesh_region> &lmr = registred_mesh_regions[&m];
- for (const mesh_region &rg : lmr)
- if (rg.compare(m, region, m)) return rg;
- lmr.push_back(region);
- return lmr.back();
- }
-
-
- void ga_workspace::add_tree(ga_tree &tree, const mesh &m,
- const mesh_im &mim, const mesh_region &rg,
- const std::string &expr,
- size_type add_derivative_order,
- bool function_expr, size_type for_interpolation,
- const std::string varname_interpolation) {
- if (tree.root) {
-
- // Eliminate the term if it corresponds to disabled variables
- if ((tree.root->test_function_type >= 1 &&
- is_disabled_variable(tree.root->name_test1)) ||
- (tree.root->test_function_type >= 2 &&
- is_disabled_variable(tree.root->name_test2))) {
- // cout<<"disabling term "; ga_print_node(tree.root, cout);
cout<<endl;
- return;
- }
- // cout << "add tree with tests functions of " << tree.root->name_test1
- // << " and " << tree.root->name_test2 << endl;
- // ga_print_node(tree.root, cout); cout << endl;
- bool remain = true;
- size_type order = 0, ind_tree = 0;
-
- if (for_interpolation)
- order = size_type(-1) - add_derivative_order;
- else {
- switch(tree.root->test_function_type) {
- case 0: order = 0; break;
- case 1: order = 1; break;
- case 3: order = 2; break;
- default: GMM_ASSERT1(false, "Inconsistent term "
- << tree.root->test_function_type);
- }
- }
-
- bool found = false;
- for (size_type i = 0; i < trees.size(); ++i) {
- if (trees[i].mim == &mim && trees[i].m == &m &&
- trees[i].order == order &&
- trees[i].name_test1.compare(tree.root->name_test1) == 0 &&
- trees[i].interpolate_name_test1.compare
- (tree.root->interpolate_name_test1) == 0 &&
- trees[i].name_test2.compare(tree.root->name_test2) == 0 &&
- trees[i].interpolate_name_test2.compare
- (tree.root->interpolate_name_test2) == 0 &&
- trees[i].rg == &rg && trees[i].interpolation == for_interpolation
&&
- trees[i].varname_interpolation.compare(varname_interpolation)==0) {
- ga_tree &ftree = *(trees[i].ptree);
-
- ftree.insert_node(ftree.root, GA_NODE_OP);
- ftree.root->op_type = GA_PLUS;
- ftree.root->children.resize(2, nullptr);
- ftree.copy_node(tree.root, ftree.root, ftree.root->children[1]);
- ga_semantic_analysis("", ftree, *this, m.dim(),
- ref_elt_dim_of_mesh(m), false, function_expr);
- found = true;
- break;
- }
- }
-
- if (!found) {
- ind_tree = trees.size(); remain = false;
- trees.push_back(tree_description());
- trees.back().mim = &mim; trees.back().m = &m;
- trees.back().rg = &rg;
- trees.back().ptree = new ga_tree;
- trees.back().ptree->swap(tree);
- pga_tree_node root = trees.back().ptree->root;
- trees.back().name_test1 = root->name_test1;
- trees.back().name_test2 = root->name_test2;
- trees.back().interpolate_name_test1 = root->interpolate_name_test1;
- trees.back().interpolate_name_test2 = root->interpolate_name_test2;
- trees.back().order = order;
- trees.back().interpolation = for_interpolation;
- trees.back().varname_interpolation = varname_interpolation;
- }
-
- if (for_interpolation == 0 && order < add_derivative_order) {
- std::set<var_trans_pair> expr_variables;
- ga_extract_variables((remain ? tree : *(trees[ind_tree].ptree)).root,
- *this, m, expr_variables, true);
- for (const var_trans_pair &var : expr_variables) {
- if (!(is_constant(var.varname))) {
- ga_tree dtree = (remain ? tree : *(trees[ind_tree].ptree));
- // cout << "Derivation with respect to " << var.first << " : "
- // << var.second << " of " << ga_tree_to_string(dtree) << endl;
- GA_TIC;
- ga_derivative(dtree, *this, m, var.varname, var.transname,
1+order);
- // cout << "Result : " << ga_tree_to_string(dtree) << endl;
- GA_TOCTIC("Derivative time");
- ga_semantic_analysis(expr, dtree, *this, m.dim(),
- ref_elt_dim_of_mesh(m), false, function_expr);
- GA_TOCTIC("Analysis after Derivative time");
- // cout << "after analysis " << ga_tree_to_string(dtree) << endl;
- add_tree(dtree, m, mim, rg, expr, add_derivative_order,
- function_expr, for_interpolation, varname_interpolation);
- }
- }
- }
- }
- }
-
- ga_workspace::m_tree::~m_tree() { if (ptree) delete ptree; }
- ga_workspace::m_tree::m_tree(const m_tree& o)
- : ptree(o.ptree), meshdim(o.meshdim), ignore_X(o.ignore_X)
- { if (o.ptree) ptree = new ga_tree(*(o.ptree)); }
- ga_workspace::m_tree &ga_workspace::m_tree::operator =(const m_tree& o) {
- if (ptree) delete ptree;
- ptree = o.ptree; meshdim = o.meshdim; ignore_X = o.ignore_X;
- if (o.ptree) ptree = new ga_tree(*(o.ptree));
- return *this;
- }
-
- size_type ga_workspace::add_expression(const std::string &expr,
- const mesh_im &mim,
- const mesh_region &rg_,
- size_type add_derivative_order) {
- const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
- // cout << "adding expression " << expr << endl;
- GA_TIC;
- size_type max_order = 0;
- std::vector<ga_tree> ltrees(1);
- ga_read_string(expr, ltrees[0]);
- // cout << "read : " << ga_tree_to_string(ltrees[0]) << endl;
- ga_semantic_analysis(expr, ltrees[0], *this, mim.linked_mesh().dim(),
- ref_elt_dim_of_mesh(mim.linked_mesh()),
- false, false, 1);
- // cout << "analysed : " << ga_tree_to_string(ltrees[0]) << endl;
- GA_TOC("First analysis time");
- if (ltrees[0].root) {
- if (test1.size() > 1 || test2.size() > 1) {
- size_type ntest2 = std::max(size_type(1), test2.size());
- size_type nb_ltrees = test1.size()*ntest2;
- ltrees.resize(nb_ltrees);
- for (size_type i = 1; i < nb_ltrees; ++i) ltrees[i] = ltrees[0];
- std::set<var_trans_pair>::iterator it1 = test1.begin();
- for (size_type i = 0; i < test1.size(); ++i, ++it1) {
- std::set<var_trans_pair>::iterator it2 = test2.begin();
- for (size_type j = 0; j < ntest2; ++j) {
- selected_test1 = *it1;
- if (test2.size()) selected_test2 = *it2++;
- // cout << "analysis with " << selected_test1.first << endl;
- ga_semantic_analysis(expr, ltrees[i*ntest2+j], *this,
- mim.linked_mesh().dim(),
- ref_elt_dim_of_mesh(mim.linked_mesh()),
- false, false, 2);
- // cout <<"split: "<< ga_tree_to_string(ltrees[i*ntest2+j]) <<
endl;
- }
- }
- }
-
- for (size_type i = 0; i < ltrees.size(); ++i) {
- if (ltrees[i].root) {
- // cout << "adding tree " << ga_tree_to_string(ltrees[i]) << endl;
- max_order = std::max(ltrees[i].root->nb_test_functions(), max_order);
- add_tree(ltrees[i], mim.linked_mesh(), mim, rg, expr,
- add_derivative_order, true, 0, "");
- }
- }
- }
- GA_TOC("Time for add expression");
- return max_order;
- }
-
- void ga_workspace::add_function_expression(const std::string &expr) {
- ga_tree tree;
- ga_read_string(expr, tree);
- ga_semantic_analysis(expr, tree, *this, 1, 1, false, true);
- if (tree.root) {
- // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
- // "Invalid function expression");
- add_tree(tree, dummy_mesh(), dummy_mesh_im(), dummy_mesh_region(),
- expr, 0, true, 0, "");
- }
- }
-
- void ga_workspace::add_interpolation_expression(const std::string &expr,
- const mesh &m,
- const mesh_region &rg_) {
- const mesh_region &rg = register_region(m, rg_);
- ga_tree tree;
- ga_read_string(expr, tree);
- ga_semantic_analysis(expr, tree, *this, m.dim(), ref_elt_dim_of_mesh(m),
- false, false);
- if (tree.root) {
- // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
- // "Invalid expression containing test functions");
- add_tree(tree, m, dummy_mesh_im(), rg, expr, 0, false, 1, "");
- }
- }
-
- void ga_workspace::add_interpolation_expression(const std::string &expr,
- const mesh_im &mim,
- const mesh_region &rg_) {
- const mesh &m = mim.linked_mesh();
- const mesh_region &rg = register_region(m, rg_);
- ga_tree tree;
- ga_read_string(expr, tree);
- ga_semantic_analysis(expr, tree, *this, m.dim(), ref_elt_dim_of_mesh(m),
- false, false);
- if (tree.root) {
- GMM_ASSERT1(tree.root->nb_test_functions() == 0,
- "Invalid expression containing test functions");
- add_tree(tree, m, mim, rg, expr, 0, false, 1, "");
- }
- }
-
- void ga_workspace::add_assignment_expression
- (const std::string &varname, const std::string &expr, const mesh_region &rg_,
- size_type order, bool before) {
- const im_data *imd = associated_im_data(varname);
- GMM_ASSERT1(imd != 0, "Only applicable to im_data");
- const mesh_im &mim = imd->linked_mesh_im();
- const mesh &m = mim.linked_mesh();
- const mesh_region &rg = register_region(m, rg_);
- ga_tree tree;
- ga_read_string(expr, tree);
- ga_semantic_analysis(expr, tree, *this, m.dim(), ref_elt_dim_of_mesh(m),
- false, false);
- if (tree.root) {
- GMM_ASSERT1(tree.root->nb_test_functions() == 0,
- "Invalid expression containing test functions");
- add_tree(tree, m, mim, rg, expr, order+1, false, (before ? 1 : 2),
- varname);
- }
- }
-
- size_type ga_workspace::nb_trees() const { return trees.size(); }
-
- ga_workspace::tree_description &ga_workspace::tree_info(size_type i)
- { return trees[i]; }
-
- bool ga_workspace::used_variables(std::vector<std::string> &vl,
- std::vector<std::string> &vl_test1,
- std::vector<std::string> &vl_test2,
- std::vector<std::string> &dl,
- size_type order) {
- bool islin = true;
- std::set<var_trans_pair> vll, dll;
- for (size_type i = 0; i < vl.size(); ++i)
- vll.insert(var_trans_pair(vl[i], ""));
- for (size_type i = 0; i < dl.size(); ++i)
- dll.insert(var_trans_pair(dl[i], ""));
-
- for (size_type i = 0; i < trees.size(); ++i) {
- ga_workspace::tree_description &td = trees[i];
- std::set<var_trans_pair> dllaux;
- bool fv = ga_extract_variables(td.ptree->root, *this, *(td.m),
- dllaux, false);
-
- if (td.order == order) {
- for (std::set<var_trans_pair>::iterator it = dllaux.begin();
- it!=dllaux.end(); ++it)
- dll.insert(*it);
- }
- switch (td.order) {
- case 0: break;
- case 1:
- if (td.order == order) {
- if (variable_group_exists(td.name_test1)) {
- for (const std::string &t : variable_group(td.name_test1))
- vll.insert(var_trans_pair(t, td.interpolate_name_test1));
- } else {
- vll.insert(var_trans_pair(td.name_test1,
- td.interpolate_name_test1));
- }
- bool found = false;
- for (const std::string &t : vl_test1)
- if (td.name_test1.compare(t) == 0)
- found = true;
- if (!found)
- vl_test1.push_back(td.name_test1);
- }
- break;
- case 2:
- if (td.order == order) {
- if (variable_group_exists(td.name_test1)) {
- for (const std::string &t : variable_group(td.name_test1))
- vll.insert(var_trans_pair(t, td.interpolate_name_test1));
- } else {
- vll.insert(var_trans_pair(td.name_test1,
- td.interpolate_name_test1));
- }
- if (variable_group_exists(td.name_test2)) {
- for (const std::string &t : variable_group(td.name_test2))
- vll.insert(var_trans_pair(t, td.interpolate_name_test2));
- } else {
- vll.insert(var_trans_pair(td.name_test2,
- td.interpolate_name_test2));
- }
- bool found = false;
- for (size_type j = 0; j < vl_test1.size(); ++j)
- if ((td.name_test1.compare(vl_test1[j]) == 0) &&
- (td.name_test2.compare(vl_test2[j]) == 0))
- found = true;
- if (!found) {
- vl_test1.push_back(td.name_test1);
- vl_test2.push_back(td.name_test2);
- }
- }
- if (fv) islin = false;
- break;
- }
- }
- vl.clear();
- for (const auto &var : vll)
- if (vl.size() == 0 || var.varname.compare(vl.back()))
- vl.push_back(var.varname);
- dl.clear();
- for (const auto &var : dll)
- if (dl.size() == 0 || var.varname.compare(dl.back()))
- dl.push_back(var.varname);
-
- return islin;
- }
-
- void ga_workspace::define_variable_group(const std::string &group_name,
- const std::vector<std::string> &nl)
{
- GMM_ASSERT1(!(variable_exists(group_name)), "The name of a group of "
- "variables cannot be the same as a variable name");
-
- std::set<const mesh *> ms;
- bool is_data_ = false;
- for (size_type i = 0; i < nl.size(); ++i) {
- if (i == 0)
- is_data_ = is_constant(nl[i]);
- else {
- GMM_ASSERT1(is_data_ == is_constant(nl[i]),
- "It is not possible to mix variables and data in a group");
- }
- GMM_ASSERT1(variable_exists(nl[i]),
- "All variables in a group have to exist in the model");
- const mesh_fem *mf = associated_mf(nl[i]);
- GMM_ASSERT1(mf, "Variables in a group should be fem variables");
- GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
- "Two variables in a group cannot share the same mesh");
- ms.insert(&(mf->linked_mesh()));
- }
- variable_groups[group_name] = nl;
- }
-
-
- const std::string &ga_workspace::variable_in_group
- (const std::string &group_name, const mesh &m) const {
- if (variable_group_exists(group_name)) {
- for (const std::string &t : variable_group(group_name))
- if (&(associated_mf(t)->linked_mesh()) == &m)
- return t;
- GMM_ASSERT1(false, "No variable in this group for the given mesh");
- } else
- return group_name;
- }
-
-
- void ga_workspace::assembly(size_type order) {
- size_type ndof;
- const ga_workspace *w = this;
- while (w->parent_workspace) w = w->parent_workspace;
- if (w->md) ndof = w->md->nb_dof(); // To eventually call actualize_sizes()
-
- GA_TIC;
- ga_instruction_set gis;
- ga_compile(*this, gis, order);
- ndof = gis.nb_dof;
- size_type max_dof = gis.max_dof;
- GA_TOCTIC("Compile time");
-
- if (order == 2) {
- if (K.use_count()) {
- gmm::clear(*K);
- gmm::resize(*K, max_dof, max_dof);
- }
- gmm::clear(unreduced_K);
- gmm::resize(unreduced_K, ndof, ndof);
- }
- if (order == 1) {
- if (V.use_count()) {
- gmm::clear(*V);
- gmm::resize(*V, max_dof);
- }
- gmm::clear(unreduced_V);
- gmm::resize(unreduced_V, ndof);
- }
- E = 0;
- GA_TOCTIC("Init time");
-
- ga_exec(gis, *this);
- GA_TOCTIC("Exec time");
-
- if (order == 1) {
- MPI_SUM_VECTOR(assembled_vector());
- MPI_SUM_VECTOR(unreduced_V);
- } else if (order == 0) {
- assembled_potential() = MPI_SUM_SCALAR(assembled_potential());
- }
-
- // Deal with reduced fems.
- if (order > 0) {
- std::set<std::string> vars_vec_done;
- std::set<std::pair<std::string, std::string> > vars_mat_done;
- for (ga_tree &tree : gis.trees) {
- if (tree.root) {
- if (order == 1) {
- const std::string &name = tree.root->name_test1;
- const std::vector<std::string> vnames_(1,name);
- const std::vector<std::string> &vnames
- = variable_group_exists(name) ? variable_group(name)
- : vnames_;
- for (const std::string &vname : vnames) {
- const mesh_fem *mf = associated_mf(vname);
- if (mf && mf->is_reduced() &&
- vars_vec_done.find(vname) == vars_vec_done.end()) {
- gmm::mult_add(gmm::transposed(mf->extension_matrix()),
- gmm::sub_vector(unreduced_V,
- gis.var_intervals[vname]),
- gmm::sub_vector(*V,
- interval_of_variable(vname)));
- vars_vec_done.insert(vname);
- }
- }
- } else {
- std::string &name1 = tree.root->name_test1;
- std::string &name2 = tree.root->name_test2;
- const std::vector<std::string> vnames1_(1,name1),
- vnames2_(2,name2);
- const std::vector<std::string> &vnames1
- = variable_group_exists(name1) ? variable_group(name1)
- : vnames1_;
- const std::vector<std::string> &vnames2
- = variable_group_exists(name2) ? variable_group(name2)
- : vnames2_;
- for (const std::string &vname1 : vnames1) {
- for (const std::string &vname2 : vnames2) {
- const mesh_fem *mf1 = associated_mf(vname1);
- const mesh_fem *mf2 = associated_mf(vname2);
- if (((mf1 && mf1->is_reduced())
- || (mf2 && mf2->is_reduced()))) {
- std::pair<std::string, std::string> p(vname1, vname2);
- if (vars_mat_done.find(p) == vars_mat_done.end()) {
- gmm::sub_interval uI1 = gis.var_intervals[vname1];
- gmm::sub_interval uI2 = gis.var_intervals[vname2];
- gmm::sub_interval I1 = interval_of_variable(vname1);
- gmm::sub_interval I2 = interval_of_variable(vname2);
- if ((mf1 && mf1->is_reduced()) &&
- (mf2 && mf2->is_reduced())) {
- model_real_sparse_matrix aux(I1.size(), uI2.size());
- model_real_row_sparse_matrix M(I1.size(), I2.size());
- gmm::mult(gmm::transposed(mf1->extension_matrix()),
- gmm::sub_matrix(unreduced_K, uI1, uI2), aux);
- gmm::mult(aux, mf2->extension_matrix(), M);
- gmm::add(M, gmm::sub_matrix(*K, I1, I2));
- } else if (mf1 && mf1->is_reduced()) {
- model_real_sparse_matrix M(I1.size(), I2.size());
- gmm::mult(gmm::transposed(mf1->extension_matrix()),
- gmm::sub_matrix(unreduced_K, uI1, uI2), M);
- gmm::add(M, gmm::sub_matrix(*K, I1, I2));
- } else {
- model_real_row_sparse_matrix M(I1.size(), I2.size());
- gmm::mult(gmm::sub_matrix(unreduced_K, uI1, uI2),
- mf2->extension_matrix(), M);
- gmm::add(M, gmm::sub_matrix(*K, I1, I2));
- }
- vars_mat_done.insert(p);
- }
- }
- }
- }
- }
- }
- }
- }
- }
-
- void ga_workspace::clear_expressions() {
- trees.clear();
- macro_trees.clear();
- }
-
- void ga_workspace::print(std::ostream &str) {
- for (size_type i = 0; i < trees.size(); ++i)
- if (trees[i].ptree->root) {
- cout << "Expression tree " << i << " of order " <<
- trees[i].ptree->root->nb_test_functions() << " :" << endl;
- ga_print_node(trees[i].ptree->root, str);
- cout << endl;
- }
- }
-
- void ga_workspace::tree_description::copy(const tree_description& td) {
- order = td.order;
- interpolation = td.interpolation;
- varname_interpolation = td.varname_interpolation;
- name_test1 = td.name_test1;
- name_test2 = td.name_test2;
- interpolate_name_test1 = td.interpolate_name_test1;
- interpolate_name_test2 = td.interpolate_name_test2;
- mim = td.mim;
- m = td.m;
- rg = td.rg;
- ptree = 0;
- if (td.ptree) ptree = new ga_tree(*(td.ptree));
- }
-
- ga_workspace::tree_description &ga_workspace::tree_description::operator =
- (const ga_workspace::tree_description& td)
- { if (ptree) delete ptree; ptree = 0; copy(td); return *this; }
- ga_workspace::tree_description::~tree_description()
- { if (ptree) delete ptree; ptree = 0; }
-
-
-
- //=========================================================================
- // Extract the constant term of degree 1 expressions
- //=========================================================================
-
- std::string ga_workspace::extract_constant_term(const mesh &m) {
- std::string constant_term;
- for (size_type i = 0; i < trees.size(); ++i) {
- ga_workspace::tree_description &td = trees[i];
-
- if (td.order == 1) {
- ga_tree local_tree = *(td.ptree);
- if (local_tree.root)
- ga_node_extract_constant_term(local_tree, local_tree.root, *this, m);
- if (local_tree.root)
- ga_semantic_analysis("", local_tree, *this, m.dim(),
- ref_elt_dim_of_mesh(m), false, false);
- if (local_tree.root && local_tree.root->node_type != GA_NODE_ZERO) {
- constant_term += "-("+ga_tree_to_string(local_tree)+")";
- }
- }
- }
- return constant_term;
- }
-
- //=========================================================================
- // Extract the order zero term
- //=========================================================================
-
- std::string ga_workspace::extract_order0_term() {
- std::string term;
- for (size_type i = 0; i < trees.size(); ++i) {
- ga_workspace::tree_description &td = trees[i];
- if (td.order == 0) {
- ga_tree &local_tree = *(td.ptree);
- if (term.size())
- term += "+("+ga_tree_to_string(local_tree)+")";
- else
- term = "("+ga_tree_to_string(local_tree)+")";
- }
- }
- return term;
- }
-
-
- //=========================================================================
- // Extract the order one term corresponding to a certain test function
- //=========================================================================
-
- std::string ga_workspace::extract_order1_term(const std::string &varname) {
- std::string term;
- for (size_type i = 0; i < trees.size(); ++i) {
- ga_workspace::tree_description &td = trees[i];
- if (td.order == 1 && td.name_test1.compare(varname) == 0) {
- ga_tree &local_tree = *(td.ptree);
- if (term.size())
- term += "+("+ga_tree_to_string(local_tree)+")";
- else
- term = "("+ga_tree_to_string(local_tree)+")";
- }
- }
- return term;
- }
-
- //=========================================================================
- // Extract Neumann terms
- //=========================================================================
-
- std::string ga_workspace::extract_Neumann_term(const std::string &varname) {
- std::string result;
- for (size_type i = 0; i < trees.size(); ++i) {
- ga_workspace::tree_description &td = trees[i];
- if (td.order == 1 && td.name_test1.compare(varname) == 0) {
- ga_tree &local_tree = *(td.ptree);
- if (local_tree.root)
- ga_extract_Neumann_term(local_tree, varname, *this,
- local_tree.root, result);
- }
- }
- return result;
- }
-
//=========================================================================
// Compilation of assembly trees into a list of basic instructions
//=========================================================================
@@ -6894,8 +5934,8 @@ namespace getfem {
}
}
- static void ga_compile_interpolation(ga_workspace &workspace,
- ga_instruction_set &gis) {
+ void ga_compile_interpolation(ga_workspace &workspace,
+ ga_instruction_set &gis) {
gis.transformations.clear();
gis.whole_instructions.clear();
for (size_type i = 0; i < workspace.nb_trees(); ++i) {
@@ -7151,9 +6191,9 @@ namespace getfem {
}
}
- static void ga_interpolation_exec(ga_instruction_set &gis,
- ga_workspace &workspace,
- ga_interpolation_context &gic) {
+ void ga_interpolation_exec(ga_instruction_set &gis,
+ ga_workspace &workspace,
+ ga_interpolation_context &gic) {
base_matrix G;
base_small_vector un, up;
@@ -7243,7 +6283,7 @@ namespace getfem {
gic.finalize();
}
- static void ga_interpolation_single_point_exec
+ void ga_interpolation_single_point_exec
(ga_instruction_set &gis, ga_workspace &workspace,
const fem_interpolation_context &ctx_x, const base_small_vector &Normal,
const mesh &interp_mesh) {
@@ -7386,954 +6426,5 @@ namespace getfem {
workspace.interpolate_transformation(t)->finalize();
}
- //=========================================================================
- // User defined functions
- //=========================================================================
-
- ga_function::ga_function(const ga_workspace &workspace_,
- const std::string &e)
- : local_workspace(true, workspace_), expr(e), gis(0) {}
-
- ga_function::ga_function(const model &md, const std::string &e)
- : local_workspace(md), expr(e), gis(0) {}
-
- ga_function::ga_function(const std::string &e)
- : local_workspace(), expr(e), gis(0) {}
-
- ga_function::ga_function(const ga_function &gaf)
- : local_workspace(gaf.local_workspace), expr(gaf.expr), gis(0)
- { if (gaf.gis) compile(); }
-
- void ga_function::compile() const {
- if (gis) delete gis;
- gis = new ga_instruction_set;
- local_workspace.clear_expressions();
- local_workspace.add_function_expression(expr);
- ga_compile_function(local_workspace, *gis, false);
- }
-
- ga_function &ga_function::operator =(const ga_function &gaf) {
- if (gis) delete gis;
- gis = 0;
- local_workspace = gaf.local_workspace;
- expr = gaf.expr;
- if (gaf.gis) compile();
- return *this;
- }
-
- ga_function::~ga_function() { if (gis) delete gis; gis = 0; }
-
- const base_tensor &ga_function::eval() const {
- GMM_ASSERT1(gis, "Uncompiled function");
- gmm::clear(local_workspace.assembled_tensor().as_vector());
- ga_function_exec(*gis);
- return local_workspace.assembled_tensor();
- }
-
- void ga_function::derivative(const std::string &var) {
- GMM_ASSERT1(gis, "Uncompiled function");
- if (local_workspace.nb_trees()) {
- ga_tree tree = *(local_workspace.tree_info(0).ptree);
- ga_derivative(tree, local_workspace, *((const mesh *)(0)), var, "", 1);
- if (tree.root) {
- ga_semantic_analysis(expr, tree, local_workspace, 1, 1, false, true);
- // To be improved to suppress test functions in the expression ...
- // ga_replace_test_by_cte do not work in all operations like
- // vector components x(1)
- // ga_replace_test_by_cte(tree.root, false);
- // ga_semantic_analysis(expr, tree, local_workspace, 1, 1,
- // false, true);
- }
- expr = ga_tree_to_string(tree);
- }
- if (gis) delete gis;
- gis = 0;
- compile();
- }
-
- //=========================================================================
- // Interpolation functions
- //=========================================================================
-
- // general Interpolation
- void ga_interpolation(ga_workspace &workspace,
- ga_interpolation_context &gic) {
- ga_instruction_set gis;
- ga_compile_interpolation(workspace, gis);
- ga_interpolation_exec(gis, workspace, gic);
- }
-
- // Interpolation on a Lagrange fem on the same mesh
- struct ga_interpolation_context_fem_same_mesh
- : public ga_interpolation_context {
- base_vector &result;
- std::vector<int> dof_count;
- const mesh_fem &mf;
- bool initialized;
- bool is_torus;
- size_type s;
-
- virtual bgeot::pstored_point_tab
- ppoints_for_element(size_type cv, short_type f,
- std::vector<size_type> &ind) const {
- pfem pf = mf.fem_of_element(cv);
- GMM_ASSERT1(pf->is_lagrange(),
- "Only Lagrange fems can be used in interpolation");
-
- if (f != short_type(-1)) {
-
- for (size_type i = 0;
- i < pf->node_convex(cv).structure()->nb_points_of_face(f); ++i)
- ind.push_back
- (pf->node_convex(cv).structure()->ind_points_of_face(f)[i]);
- } else {
- for (size_type i = 0; i < pf->node_convex(cv).nb_points(); ++i)
- ind.push_back(i);
- }
-
- return pf->node_tab(cv);
- }
-
- virtual bool use_pgp(size_type) const { return true; }
- virtual bool use_mim() const { return false; }
-
- void init_(size_type si, size_type q, size_type qmult) {
- s = si;
- gmm::resize(result, qmult * mf.nb_basic_dof());
- gmm::clear(result);
- gmm::resize(dof_count, mf.nb_basic_dof()/q);
- gmm::clear(dof_count);
- initialized = true;
- }
-
- void store_result_for_torus(size_type cv, size_type i, base_tensor &t) {
- size_type target_dim = mf.fem_of_element(cv)->target_dim();
- GMM_ASSERT2(target_dim == 3, "Invalid torus fem.");
- size_type qdim = 1;
- size_type result_dim = 2;
- if (!initialized) {init_(qdim, qdim, qdim);}
- size_type idof = mf.ind_basic_dof_of_element(cv)[i];
- result[idof] = t[idof%result_dim];
- ++dof_count[idof];
- }
-
- virtual void store_result(size_type cv, size_type i, base_tensor &t) {
- if (is_torus){
- store_result_for_torus(cv, i, t);
- return;
- }
- size_type si = t.size();
- size_type q = mf.get_qdim();
- size_type qmult = si / q;
- GMM_ASSERT1( (si % q) == 0, "Incompatibility between the mesh_fem and "
- "the size of the expression to be interpolated");
- if (!initialized) { init_(si, q, qmult); }
- GMM_ASSERT1(s == si, "Internal error");
- size_type idof = mf.ind_basic_dof_of_element(cv)[i*q];
- gmm::add(t.as_vector(),
- gmm::sub_vector(result, gmm::sub_interval(qmult*idof, s)));
- (dof_count[idof/q])++;
- }
-
- virtual void finalize() {
- std::vector<size_type> data(3);
- data[0] = initialized ? result.size() : 0;
- data[1] = initialized ? dof_count.size() : 0;
- data[2] = initialized ? s : 0;
- MPI_MAX_VECTOR(data);
- if (!initialized) {
- if (data[0]) {
- gmm::resize(result, data[0]);
- gmm::resize(dof_count, data[1]);
- gmm::clear(dof_count);
- s = data[2];
- }
- gmm::clear(result);
- }
- if (initialized)
- GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[2] &&
- gmm::vect_size(dof_count) == data[1], "Incompatible sizes");
- MPI_SUM_VECTOR(result);
- MPI_SUM_VECTOR(dof_count);
- for (size_type i = 0; i < dof_count.size(); ++i)
- if (dof_count[i])
- gmm::scale(gmm::sub_vector(result, gmm::sub_interval(s*i, s)),
- scalar_type(1) / scalar_type(dof_count[i]));
- }
-
- virtual const mesh &linked_mesh() { return mf.linked_mesh(); }
-
- ga_interpolation_context_fem_same_mesh(const mesh_fem &mf_, base_vector &r)
- : result(r), mf(mf_), initialized(false), is_torus(false) {
- GMM_ASSERT1(!(mf.is_reduced()),
- "Interpolation on reduced fem is not allowed");
- if (dynamic_cast<const torus_mesh_fem*>(&mf)){
- auto first_cv = mf.first_convex_of_basic_dof(0);
- auto target_dim = mf.fem_of_element(first_cv)->target_dim();
- if (target_dim == 3) is_torus = true;
- }
- }
- };
-
- void ga_interpolation_Lagrange_fem
- (ga_workspace &workspace, const mesh_fem &mf, base_vector &result) {
- ga_interpolation_context_fem_same_mesh gic(mf, result);
- ga_interpolation(workspace, gic);
- }
-
- void ga_interpolation_Lagrange_fem
- (const getfem::model &md, const std::string &expr, const mesh_fem &mf,
- base_vector &result, const mesh_region &rg) {
- ga_workspace workspace(md);
- workspace.add_interpolation_expression(expr, mf.linked_mesh(), rg);
- ga_interpolation_Lagrange_fem(workspace, mf, result);
- }
-
- // Interpolation on a cloud of points
- struct ga_interpolation_context_mti
- : public ga_interpolation_context {
- base_vector &result;
- const mesh_trans_inv &mti;
- bool initialized;
- size_type s, nbdof;
-
-
- virtual bgeot::pstored_point_tab
- ppoints_for_element(size_type cv, short_type,
- std::vector<size_type> &ind) const {
- std::vector<size_type> itab;
- mti.points_on_convex(cv, itab);
- std::vector<base_node> pt_tab(itab.size());
- for (size_type i = 0; i < itab.size(); ++i) {
- pt_tab[i] = mti.reference_coords()[itab[i]];
- ind.push_back(i);
- }
- return store_point_tab(pt_tab);
- }
-
- virtual bool use_pgp(size_type) const { return false; }
- virtual bool use_mim() const { return false; }
-
- virtual void store_result(size_type cv, size_type i, base_tensor &t) {
- size_type si = t.size();
- if (!initialized) {
- s = si;
- gmm::resize(result, s * nbdof);
- gmm::clear(result);
- initialized = true;
- }
- GMM_ASSERT1(s == si, "Internal error");
- size_type ipt = mti.point_on_convex(cv, i);
- size_type dof_t = mti.id_of_point(ipt);
- gmm::add(t.as_vector(),
- gmm::sub_vector(result, gmm::sub_interval(s*dof_t, s)));
- }
-
- virtual void finalize() {
- std::vector<size_type> data(2);
- data[0] = initialized ? result.size() : 0;
- data[1] = initialized ? s : 0;
- MPI_MAX_VECTOR(data);
- if (!initialized) {
- if (data[0]) {
- gmm::resize(result, data[0]);
- s = data[1];
- }
- gmm::clear(result);
- }
- if (initialized)
- GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[1],
- "Incompatible sizes");
- MPI_SUM_VECTOR(result);
- }
-
- virtual const mesh &linked_mesh() { return mti.linked_mesh(); }
-
- ga_interpolation_context_mti(const mesh_trans_inv &mti_, base_vector &r,
- size_type nbdof_ = size_type(-1))
- : result(r), mti(mti_), initialized(false), nbdof(nbdof_) {
- if (nbdof == size_type(-1)) nbdof = mti.nb_points();
- }
- };
-
- // Distribute to be parallelized
- void ga_interpolation_mti
- (const getfem::model &md, const std::string &expr, mesh_trans_inv &mti,
- base_vector &result, const mesh_region &rg, int extrapolation,
- const mesh_region &rg_source, size_type nbdof) {
-
- ga_workspace workspace(md);
- workspace.add_interpolation_expression(expr, mti.linked_mesh(), rg);
-
- mti.distribute(extrapolation, rg_source);
- ga_interpolation_context_mti gic(mti, result, nbdof);
- ga_interpolation(workspace, gic);
- }
-
-
- // Interpolation on a im_data
- struct ga_interpolation_context_im_data
- : public ga_interpolation_context {
- base_vector &result;
- const im_data &imd;
- bool initialized;
- size_type s;
-
- virtual bgeot::pstored_point_tab
- ppoints_for_element(size_type cv, short_type f,
- std::vector<size_type> &ind) const {
- pintegration_method pim =imd.linked_mesh_im().int_method_of_element(cv);
- if (pim->type() == IM_NONE) return bgeot::pstored_point_tab();
- GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot "
- "be used in high level generic assembly");
- size_type i_start(0), i_end(0);
- if (f == short_type(-1))
- i_end = pim->approx_method()->nb_points_on_convex();
- else {
- i_start = pim->approx_method()->ind_first_point_on_face(f);
- i_end = i_start + pim->approx_method()->nb_points_on_face(f);
- }
- for (size_type i = i_start; i < i_end; ++i) ind.push_back(i);
- return pim->approx_method()->pintegration_points();
- }
-
- virtual bool use_pgp(size_type cv) const {
- pintegration_method pim =imd.linked_mesh_im().int_method_of_element(cv);
- if (pim->type() == IM_NONE) return false;
- GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot "
- "be used in high level generic assembly");
- return !(pim->approx_method()->is_built_on_the_fly());
- }
- virtual bool use_mim() const { return true; }
-
- virtual void store_result(size_type cv, size_type i, base_tensor &t) {
- size_type si = t.size();
- if (!initialized) {
- s = si;
- GMM_ASSERT1(imd.tensor_size() == t.sizes() ||
- (imd.tensor_size().size() == size_type(1) &&
- imd.tensor_size()[0] == size_type(1) &&
- si == size_type(1)),
- "Im_data tensor size " << imd.tensor_size() <<
- " does not match the size of the interpolated "
- "expression " << t.sizes() << ".");
- gmm::resize(result, s * imd.nb_filtered_index());
- gmm::clear(result);
- initialized = true;
- }
- GMM_ASSERT1(s == si, "Internal error");
- size_type ipt = imd.filtered_index_of_point(cv, i);
- GMM_ASSERT1(ipt != size_type(-1),
- "Im data with no data on the current integration point.");
- gmm::add(t.as_vector(),
- gmm::sub_vector(result, gmm::sub_interval(s*ipt, s)));
- }
-
- virtual void finalize() {
- std::vector<size_type> data(2);
- data[0] = initialized ? result.size() : 0;
- data[1] = initialized ? s : 0;
- MPI_MAX_VECTOR(data);
- if (initialized) {
- GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[1],
- "Incompatible sizes");
- } else {
- if (data[0]) {
- gmm::resize(result, data[0]);
- s = data[1];
- }
- gmm::clear(result);
- }
- MPI_SUM_VECTOR(result);
- }
-
- virtual const mesh &linked_mesh() { return imd.linked_mesh(); }
-
- ga_interpolation_context_im_data(const im_data &imd_, base_vector &r)
- : result(r), imd(imd_), initialized(false) { }
- };
-
- void ga_interpolation_im_data
- (ga_workspace &workspace, const im_data &imd, base_vector &result) {
- ga_interpolation_context_im_data gic(imd, result);
- ga_interpolation(workspace, gic);
- }
-
- void ga_interpolation_im_data
- (const getfem::model &md, const std::string &expr, const im_data &imd,
- base_vector &result, const mesh_region &rg) {
- ga_workspace workspace(md);
- workspace.add_interpolation_expression
- (expr, imd.linked_mesh_im(), rg);
-
- ga_interpolation_im_data(workspace, imd, result);
- }
-
-
- // Interpolation on a stored_mesh_slice
- struct ga_interpolation_context_mesh_slice
- : public ga_interpolation_context {
- base_vector &result;
- const stored_mesh_slice &sl;
- bool initialized;
- size_type s;
- std::vector<size_type> first_node;
-
- virtual bgeot::pstored_point_tab
- ppoints_for_element(size_type cv, short_type f,
- std::vector<size_type> &ind) const {
- GMM_ASSERT1(f == short_type(-1), "No support for interpolation on faces"
- " for a stored_mesh_slice yet.");
- size_type ic = sl.convex_pos(cv);
- const mesh_slicer::cs_nodes_ct &nodes = sl.nodes(ic);
- std::vector<base_node> pt_tab(nodes.size());
- for (size_type i=0; i < nodes.size(); ++i) {
- pt_tab[i] = nodes[i].pt_ref;
- ind.push_back(i);
- }
- return store_point_tab(pt_tab);
- }
-
- virtual bool use_pgp(size_type /* cv */) const { return false; } // why
not?
- virtual bool use_mim() const { return false; }
-
- virtual void store_result(size_type cv, size_type i, base_tensor &t) {
- size_type si = t.size();
- if (!initialized) {
- s = si;
- gmm::resize(result, s * sl.nb_points());
- gmm::clear(result);
- initialized = true;
- first_node.resize(sl.nb_convex());
- for (size_type ic=0; ic < sl.nb_convex()-1; ++ic)
- first_node[ic+1] = first_node[ic] + sl.nodes(ic).size();
- }
- GMM_ASSERT1(s == si && result.size() == s * sl.nb_points(), "Internal
error");
- size_type ic = sl.convex_pos(cv);
- size_type ipt = first_node[ic] + i;
- gmm::add(t.as_vector(),
- gmm::sub_vector(result, gmm::sub_interval(s*ipt, s)));
- }
-
- virtual void finalize() {
- std::vector<size_type> data(2);
- data[0] = initialized ? result.size() : 0;
- data[1] = initialized ? s : 0;
- MPI_MAX_VECTOR(data);
- if (initialized) {
- GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[1],
- "Incompatible sizes");
- } else {
- if (data[0]) {
- gmm::resize(result, data[0]);
- s = data[1];
- }
- gmm::clear(result);
- }
- MPI_SUM_VECTOR(result);
- }
-
- virtual const mesh &linked_mesh() { return sl.linked_mesh(); }
-
- ga_interpolation_context_mesh_slice(const stored_mesh_slice &sl_,
base_vector &r)
- : result(r), sl(sl_), initialized(false) { }
- };
-
- void ga_interpolation_mesh_slice
- (ga_workspace &workspace, const stored_mesh_slice &sl, base_vector &result) {
- ga_interpolation_context_mesh_slice gic(sl, result);
- ga_interpolation(workspace, gic);
- }
-
- void ga_interpolation_mesh_slice
- (const getfem::model &md, const std::string &expr, const stored_mesh_slice
&sl,
- base_vector &result, const mesh_region &rg) {
- ga_workspace workspace(md);
- workspace.add_interpolation_expression(expr, sl.linked_mesh(), rg);
- ga_interpolation_mesh_slice(workspace, sl, result);
- }
-
-
- //=========================================================================
- // Local projection functions
- //=========================================================================
-
- void ga_local_projection(const getfem::model &md, const mesh_im &mim,
- const std::string &expr, const mesh_fem &mf,
- base_vector &result, const mesh_region ®ion) {
-
- // The case where the expression is a vector one and mf a scalar fem is
- // not taken into account for the moment.
-
- // Could be improved by not performing the assembly of the global mass
matrix
- // working locally. This means a specific assembly.
- model_real_sparse_matrix M(mf.nb_dof(), mf.nb_dof());
- asm_mass_matrix(M, mim, mf, region);
-
- ga_workspace workspace(md);
- size_type nbdof = md.nb_dof();
- gmm::sub_interval I(nbdof, mf.nb_dof());
- workspace.add_fem_variable("c__dummy_var_95_", mf, I, base_vector(nbdof));
- if (mf.get_qdims().size() > 1)
-
workspace.add_expression("("+expr+"):Test_c__dummy_var_95_",mim,region,2);
- else
-
workspace.add_expression("("+expr+").Test_c__dummy_var_95_",mim,region,2);
- base_vector residual(nbdof+mf.nb_dof());
- workspace.set_assembled_vector(residual);
- workspace.assembly(1);
- getfem::base_vector F(mf.nb_dof());
- gmm::resize(result, mf.nb_dof());
- gmm::copy(gmm::sub_vector(residual, I), F);
-
- getfem::base_matrix loc_M;
- getfem::base_vector loc_U;
- for (mr_visitor v(region, mf.linked_mesh(), true); !v.finished(); ++v) {
- if (mf.convex_index().is_in(v.cv())) {
- size_type nd = mf.nb_basic_dof_of_element(v.cv());
- loc_M.base_resize(nd, nd); gmm::resize(loc_U, nd);
- gmm::sub_index J(mf.ind_basic_dof_of_element(v.cv()));
- gmm::copy(gmm::sub_matrix(M, J, J), loc_M);
- gmm::lu_solve(loc_M, loc_U, gmm::sub_vector(F, J));
- gmm::copy(loc_U, gmm::sub_vector(result, J));
- }
- }
- MPI_SUM_VECTOR(result);
- }
-
- //=========================================================================
- // Interpolate transformation with an expression
- //=========================================================================
-
- class interpolate_transformation_expression
- : public virtual_interpolate_transformation, public context_dependencies {
-
- struct workspace_gis_pair : public std::pair<ga_workspace,
ga_instruction_set> {
- inline ga_workspace &workspace() { return this->first; }
- inline ga_instruction_set &gis() { return this->second; }
- };
-
- const mesh &source_mesh;
- const mesh &target_mesh;
- std::string expr;
- mutable bgeot::rtree element_boxes;
- mutable bool recompute_elt_boxes;
- mutable ga_workspace local_workspace;
- mutable ga_instruction_set local_gis;
- mutable bgeot::geotrans_inv_convex gic;
- mutable base_node P;
- mutable std::set<var_trans_pair> used_vars;
- mutable std::set<var_trans_pair> used_data;
- mutable std::map<var_trans_pair,
- workspace_gis_pair> compiled_derivatives;
- mutable bool extract_variable_done;
- mutable bool extract_data_done;
-
- public:
- void update_from_context() const {
- recompute_elt_boxes = true;
- }
-
- void extract_variables(const ga_workspace &workspace,
- std::set<var_trans_pair> &vars,
- bool ignore_data, const mesh &/* m */,
- const std::string &/* interpolate_name */) const {
- if ((ignore_data && !extract_variable_done) ||
- (!ignore_data && !extract_data_done)) {
- if (ignore_data)
- used_vars.clear();
- else
- used_data.clear();
- ga_workspace aux_workspace;
- aux_workspace = ga_workspace(true, workspace);
- aux_workspace.clear_expressions();
- aux_workspace.add_interpolation_expression(expr, source_mesh);
- for (size_type i = 0; i < aux_workspace.nb_trees(); ++i)
- ga_extract_variables(aux_workspace.tree_info(i).ptree->root,
- aux_workspace, source_mesh,
- ignore_data ? used_vars : used_data,
- ignore_data);
- if (ignore_data)
- extract_variable_done = true;
- else
- extract_data_done = true;
- }
- if (ignore_data)
- vars.insert(used_vars.begin(), used_vars.end());
- else
- vars.insert(used_data.begin(), used_data.end());
- }
-
- void init(const ga_workspace &workspace) const {
- size_type N = target_mesh.dim();
-
- // Expression compilation
- local_workspace = ga_workspace(true, workspace);
- local_workspace.clear_expressions();
-
- local_workspace.add_interpolation_expression(expr, source_mesh);
- local_gis = ga_instruction_set();
- ga_compile_interpolation(local_workspace, local_gis);
-
- // In fact, transformations are not allowed ... for future compatibility
- for (const std::string &transname : local_gis.transformations)
- local_workspace.interpolate_transformation(transname)
- ->init(local_workspace);
-
- if (!extract_variable_done) {
- std::set<var_trans_pair> vars;
- extract_variables(workspace, vars, true, source_mesh, "");
- }
-
- for (const var_trans_pair &var : used_vars) {
- workspace_gis_pair &pwi = compiled_derivatives[var];
- pwi.workspace() = local_workspace;
- pwi.gis() = ga_instruction_set();
- if (pwi.workspace().nb_trees()) {
- ga_tree tree = *(pwi.workspace().tree_info(0).ptree);
- ga_derivative(tree, pwi.workspace(), source_mesh,
- var.varname, var.transname, 1);
- if (tree.root)
- ga_semantic_analysis(expr, tree, local_workspace, 1, 1,
- false, true);
- ga_compile_interpolation(pwi.workspace(), pwi.gis());
- }
- }
-
- // Element_boxes update (if necessary)
- if (recompute_elt_boxes) {
-
- element_boxes.clear();
- base_node bmin(N), bmax(N);
- for (dal::bv_visitor cv(target_mesh.convex_index());
- !cv.finished(); ++cv) {
-
- bgeot::pgeometric_trans pgt = target_mesh.trans_of_convex(cv);
-
- size_type nbd_t = pgt->nb_points();
- if (nbd_t) {
- gmm::copy(target_mesh.points_of_convex(cv)[0], bmin);
- gmm::copy(bmin, bmax);
- } else {
- gmm::clear(bmin);
- gmm::clear(bmax);
- }
- for (short_type ip = 1; ip < nbd_t; ++ip) {
- // size_type ind = target_mesh.ind_points_of_convex(cv)[ip];
- const base_node &pt = target_mesh.points_of_convex(cv)[ip];
-
- for (size_type k = 0; k < N; ++k) {
- bmin[k] = std::min(bmin[k], pt[k]);
- bmax[k] = std::max(bmax[k], pt[k]);
- }
- }
-
- scalar_type h = bmax[0] - bmin[0];
- for (size_type k = 1; k < N; ++k) h = std::max(h, bmax[k]-bmin[k]);
- if (pgt->is_linear()) h *= 1E-4;
- for (auto &&val : bmin) val -= h*0.2;
- for (auto &&val : bmax) val += h*0.2;
-
- element_boxes.add_box(bmin, bmax, cv);
- }
- recompute_elt_boxes = false;
- }
- }
-
- void finalize() const {
- for (const std::string &transname : local_gis.transformations)
- local_workspace.interpolate_transformation(transname)->finalize();
- local_gis = ga_instruction_set();
- }
-
-
- int transform(const ga_workspace &/*workspace*/, const mesh &m,
- fem_interpolation_context &ctx_x,
- const base_small_vector &Normal,
- const mesh **m_t,
- size_type &cv, short_type &face_num,
- base_node &P_ref,
- base_small_vector &/*N_y*/,
- std::map<var_trans_pair, base_tensor> &derivatives,
- bool compute_derivatives) const {
- int ret_type = 0;
-
- ga_interpolation_single_point_exec(local_gis, local_workspace, ctx_x,
- Normal, m);
-
- GMM_ASSERT1(local_workspace.assembled_tensor().size() == m.dim(),
- "Wrong dimension of the tranformation expression");
- P.resize(m.dim());
- gmm::copy(local_workspace.assembled_tensor().as_vector(), P);
-
- bgeot::rtree::pbox_set bset;
- element_boxes.find_boxes_at_point(P, bset);
- *m_t = &target_mesh;
-
- while (bset.size()) {
- bgeot::rtree::pbox_set::iterator it = bset.begin(), itmax = it;
-
- if (bset.size() > 1) {
- // Searching the box for which the point is the most in the interior
- scalar_type rate_max = scalar_type(-1);
- for (; it != bset.end(); ++it) {
-
- scalar_type rate_box = scalar_type(1);
- for (size_type i = 0; i < m.dim(); ++i) {
- scalar_type h = (*it)->max[i] - (*it)->min[i];
- if (h > scalar_type(0)) {
- scalar_type rate
- = std::min((*it)->max[i] - P[i], P[i] - (*it)->min[i]) / h;
- rate_box = std::min(rate, rate_box);
- }
- }
- if (rate_box > rate_max) {
- itmax = it;
- rate_max = rate_box;
- }
- }
- }
-
- cv = (*itmax)->id;
- gic.init(target_mesh.points_of_convex(cv),
- target_mesh.trans_of_convex(cv));
-
- bool converged = true;
- bool is_in = gic.invert(P, P_ref, converged, 1E-4);
- // cout << "cv = " << cv << " P = " << P << " P_ref = " << P_ref <<
endl;
- // cout << " is_in = " << int(is_in) << endl;
- // for (size_type iii = 0;
- // iii < target_mesh.points_of_convex(cv).size(); ++iii)
- // cout << target_mesh.points_of_convex(cv)[iii] << endl;
-
- if (is_in && converged) {
- face_num = short_type(-1); // Should detect potential faces ?
- ret_type = 1;
- break;
- }
-
- if (bset.size() == 1) break;
- bset.erase(itmax);
- }
-
- // Note on derivatives of the transformation : for efficiency and
- // simplicity reasons, the derivative should be computed with
- // the value of corresponding test functions. This means that
- // for a transformation F(u) the computed derivative is F'(u).Test_u
- // including the Test_u.
- if (compute_derivatives) { // To be tested both with the computation
- // of derivative. Could be optimized ?
- for (auto &&d : derivatives) {
- workspace_gis_pair &pwi = compiled_derivatives[d.first];
-
- gmm::clear(pwi.workspace().assembled_tensor().as_vector());
- ga_function_exec(pwi.gis());
- d.second = pwi.workspace().assembled_tensor();
- }
- }
- return ret_type;
- }
-
- interpolate_transformation_expression(const mesh &sm, const mesh &tm,
- const std::string &expr_)
- : source_mesh(sm), target_mesh(tm), expr(expr_),
- recompute_elt_boxes(true), extract_variable_done(false),
- extract_data_done(false)
- { this->add_dependency(tm); }
-
- };
-
-
- void add_interpolate_transformation_from_expression
- (model &md, const std::string &name, const mesh &sm, const mesh &tm,
- const std::string &expr) {
- pinterpolate_transformation
- p = std::make_shared<interpolate_transformation_expression>(sm, tm,
expr);
- md.add_interpolate_transformation(name, p);
- }
-
- void add_interpolate_transformation_from_expression
- (ga_workspace &workspace, const std::string &name, const mesh &sm,
- const mesh &tm, const std::string &expr) {
- pinterpolate_transformation
- p = std::make_shared<interpolate_transformation_expression>(sm, tm,
expr);
- workspace.add_interpolate_transformation(name, p);
- }
-
- //=========================================================================
- // Interpolate transformation on neighbour element (for internal faces)
- //=========================================================================
-
- class interpolate_transformation_neighbour
- : public virtual_interpolate_transformation, public context_dependencies {
-
- public:
- void update_from_context() const {}
- void extract_variables(const ga_workspace &/* workspace */,
- std::set<var_trans_pair> &/* vars */,
- bool /* ignore_data */, const mesh &/* m */,
- const std::string &/* interpolate_name */) const {}
- void init(const ga_workspace &/* workspace */) const {}
- void finalize() const {}
-
- int transform(const ga_workspace &/*workspace*/, const mesh &m_x,
- fem_interpolation_context &ctx_x,
- const base_small_vector &/*Normal*/, const mesh **m_t,
- size_type &cv, short_type &face_num,
- base_node &P_ref,
- base_small_vector &/*N_y*/,
- std::map<var_trans_pair, base_tensor> &/*derivatives*/,
- bool compute_derivatives) const {
-
- int ret_type = 0;
- *m_t = &m_x;
- size_type cv_x = ctx_x.convex_num();
- short_type face_x = ctx_x.face_num();
- GMM_ASSERT1(face_x != short_type(-1), "Neighbour transformation can "
- "only be applied to internal faces");
-
- auto adj_face = m_x.adjacent_face(cv_x, face_x);
-
- if (adj_face.cv != size_type(-1)) {
- bgeot::geotrans_inv_convex gic;
- gic.init(m_x.points_of_convex(adj_face.cv),
- m_x.trans_of_convex(adj_face.cv));
- bool converged = true;
- bool is_in = gic.invert(ctx_x.xreal(), P_ref, converged, 1E-4);
- GMM_ASSERT1(is_in && converged, "Geometric transformation inversion "
- "has failed in neighbour transformation");
- face_num = adj_face.f;
- cv = adj_face.cv;
- ret_type = 1;
- }
- GMM_ASSERT1(!compute_derivatives,
- "No derivative for this transformation");
- return ret_type;
- }
-
- interpolate_transformation_neighbour() { }
-
- };
-
-
- pinterpolate_transformation interpolate_transformation_neighbour_instance()
- {
- pinterpolate_transformation
- p = std::make_shared<interpolate_transformation_neighbour>();
- return p;
- }
-
- //=========================================================================
- // Interpolate transformation on neighbour element (for extrapolation)
- //=========================================================================
-
- class interpolate_transformation_element_extrapolation
- : public virtual_interpolate_transformation, public context_dependencies {
-
- const mesh &sm;
- std::map<size_type, size_type> elt_corr;
-
- public:
- void update_from_context() const {}
- void extract_variables(const ga_workspace &/* workspace */,
- std::set<var_trans_pair> &/* vars */,
- bool /* ignore_data */, const mesh &/* m */,
- const std::string &/* interpolate_name */) const {}
- void init(const ga_workspace &/* workspace */) const {}
- void finalize() const {}
-
- int transform(const ga_workspace &/*workspace*/, const mesh &m_x,
- fem_interpolation_context &ctx_x,
- const base_small_vector &/*Normal*/, const mesh **m_t,
- size_type &cv, short_type &face_num,
- base_node &P_ref,
- base_small_vector &/*N_y*/,
- std::map<var_trans_pair, base_tensor> &/*derivatives*/,
- bool compute_derivatives) const {
- int ret_type = 0;
- *m_t = &m_x;
- GMM_ASSERT1(&sm == &m_x, "Bad mesh");
- size_type cv_x = ctx_x.convex_num(), cv_y = cv_x;
- auto it = elt_corr.find(cv_x);
- if (it != elt_corr.end()) cv_y = it->second;
-
- if (cv_x != cv_y) {
- bgeot::geotrans_inv_convex gic;
- gic.init(m_x.points_of_convex(cv_y),
- m_x.trans_of_convex(cv_y));
- bool converged = true;
- gic.invert(ctx_x.xreal(), P_ref, converged, 1E-4);
- GMM_ASSERT1(converged, "Geometric transformation inversion "
- "has failed in element extrapolation transformation");
- face_num = short_type(-1);
- cv = cv_y;
- ret_type = 1;
- } else {
- cv = cv_x;
- face_num = short_type(-1);
- P_ref = ctx_x.xref();
- ret_type = 1;
- }
- GMM_ASSERT1(!compute_derivatives,
- "No derivative for this transformation");
- return ret_type;
- }
-
- void set_correspondance(const std::map<size_type, size_type> &ec) {
- elt_corr = ec;
- }
-
- interpolate_transformation_element_extrapolation
- (const mesh &sm_, const std::map<size_type, size_type> &ec)
- : sm(sm_), elt_corr(ec) { }
- };
-
-
- void add_element_extrapolation_transformation
- (model &md, const std::string &name, const mesh &sm,
- std::map<size_type, size_type> &elt_corr) {
- pinterpolate_transformation
- p = std::make_shared<interpolate_transformation_element_extrapolation>
- (sm, elt_corr);
- md.add_interpolate_transformation(name, p);
- }
-
- void add_element_extrapolation_transformation
- (ga_workspace &workspace, const std::string &name, const mesh &sm,
- std::map<size_type, size_type> &elt_corr) {
- pinterpolate_transformation
- p = std::make_shared<interpolate_transformation_element_extrapolation>
- (sm, elt_corr);
- workspace.add_interpolate_transformation(name, p);
- }
-
- void set_element_extrapolation_correspondance
- (ga_workspace &workspace, const std::string &name,
- std::map<size_type, size_type> &elt_corr) {
- GMM_ASSERT1(workspace.interpolate_transformation_exists(name),
- "Unknown transformation");
- auto pit=workspace.interpolate_transformation(name).get();
- auto cpext
- = dynamic_cast<const interpolate_transformation_element_extrapolation *>
- (pit);
- GMM_ASSERT1(cpext,
- "The transformation is not of element extrapolation type");
- const_cast<interpolate_transformation_element_extrapolation *>(cpext)
- ->set_correspondance(elt_corr);
- }
-
- void set_element_extrapolation_correspondance
- (model &md, const std::string &name,
- std::map<size_type, size_type> &elt_corr) {
- GMM_ASSERT1(md.interpolate_transformation_exists(name),
- "Unknown transformation");
- auto pit=md.interpolate_transformation(name).get();
- auto cpext
- = dynamic_cast<const interpolate_transformation_element_extrapolation *>
- (pit);
- GMM_ASSERT1(cpext,
- "The transformation is not of element extrapolation type");
- const_cast<interpolate_transformation_element_extrapolation *>(cpext)
- ->set_correspondance(elt_corr);
- }
} /* end of namespace */
diff --git a/src/getfem_generic_assembly_functions_and_operators.cc
b/src/getfem_generic_assembly_functions_and_operators.cc
index 05f6d85..36abf9b 100644
--- a/src/getfem_generic_assembly_functions_and_operators.cc
+++ b/src/getfem_generic_assembly_functions_and_operators.cc
@@ -43,23 +43,13 @@ BoostMathFunction const erfc = boost::math::erfc<double>;
namespace getfem {
- extern bool predef_operators_nonlinear_elasticity_initialized;
- extern bool predef_operators_plasticity_initialized;
- extern bool predef_operators_contact_initialized;
-
base_matrix& __mat_aux1()
{
DEFINE_STATIC_THREAD_LOCAL(base_matrix, m);
return m;
}
- //=========================================================================
- // Structure dealing with predefined special functions
- // such as mesh_dim, pi, qdim ...
- //=========================================================================
- typedef std::set<std::string> ga_spec_function_tab;
- static ga_spec_function_tab SPEC_FUNCTIONS;
//=========================================================================
// Structure dealing with predefined scalar functions.
@@ -402,12 +392,28 @@ namespace getfem {
// Initialization of predefined functions and operators.
//=========================================================================
+ ga_predef_function::ga_predef_function()
+ : expr_(""), derivative1_(""), derivative2_(""), gis(nullptr) {}
+
+ ga_predef_function::ga_predef_function(pscalar_func_onearg f,
+ size_type dtype__,
+ const std::string &der)
+ : ftype_(0), dtype_(dtype__), nbargs_(1), f1_(f), expr_(""),
+ derivative1_(der), derivative2_("") {}
+ ga_predef_function::ga_predef_function(pscalar_func_twoargs f,
+ size_type dtype__,
+ const std::string &der1,
+ const std::string &der2)
+ : ftype_(0), dtype_(dtype__), nbargs_(2), f2_(f),
+ expr_(""), derivative1_(der1), derivative2_(der2), gis(nullptr) {}
+ ga_predef_function::ga_predef_function(const std::string &expr__)
+ : ftype_(1), dtype_(3), nbargs_(1), expr_(expr__),
+ derivative1_(""), derivative2_(""), t(1, 0.), u(1, 0.), gis(nullptr) {}
- bool init_predef_functions() {
-
- // Predefined functions
- ga_predef_function_tab &PREDEF_FUNCTIONS
- = dal::singleton<ga_predef_function_tab>::instance(0);
+
+ ga_predef_function_tab::ga_predef_function_tab() {
+
+ ga_predef_function_tab &PREDEF_FUNCTIONS = *this;
// Power functions and their derivatives
PREDEF_FUNCTIONS["sqrt"] = ga_predef_function(sqrt, 1, "DER_PDFUNC_SQRT");
@@ -533,29 +539,33 @@ namespace getfem {
PREDEF_FUNCTIONS["DER_PDFUNC1_MAX"] = ga_predef_function(ga_der_t_max);
PREDEF_FUNCTIONS["DER_PDFUNC2_MAX"] = ga_predef_function(ga_der_u_max);
+ }
+ ga_spec_function_tab::ga_spec_function_tab() {
// Predefined special functions
-
+ ga_spec_function_tab &SPEC_FUNCTIONS = *this;
+
SPEC_FUNCTIONS.insert("pi");
SPEC_FUNCTIONS.insert("meshdim");
SPEC_FUNCTIONS.insert("timestep");
SPEC_FUNCTIONS.insert("qdim");
SPEC_FUNCTIONS.insert("qdims");
SPEC_FUNCTIONS.insert("Id");
+ }
+
+ ga_predef_operator_tab::ga_predef_operator_tab(void) {
// Predefined operators
- ga_predef_operator_tab &PREDEF_OPERATORS
- = dal::singleton<ga_predef_operator_tab>::instance();
+ ga_predef_operator_tab &PREDEF_OPERATORS = *this;
PREDEF_OPERATORS.add_method("Norm", std::make_shared<norm_operator>());
PREDEF_OPERATORS.add_method("Norm_sqr",
std::make_shared<norm_sqr_operator>());
PREDEF_OPERATORS.add_method("Det", std::make_shared<det_operator>());
PREDEF_OPERATORS.add_method("Inv", std::make_shared<inverse_operator>());
- return true;
}
- bool predef_functions_initialized = init_predef_functions();
+
bool ga_function_exists(const std::string &name) {
const ga_predef_function_tab &PREDEF_FUNCTIONS
@@ -642,4 +652,69 @@ namespace getfem {
}
}
+ //=========================================================================
+ // User defined functions
+ //=========================================================================
+
+ ga_function::ga_function(const ga_workspace &workspace_,
+ const std::string &e)
+ : local_workspace(true, workspace_), expr(e), gis(0) {}
+
+ ga_function::ga_function(const model &md, const std::string &e)
+ : local_workspace(md), expr(e), gis(0) {}
+
+ ga_function::ga_function(const std::string &e)
+ : local_workspace(), expr(e), gis(0) {}
+
+ ga_function::ga_function(const ga_function &gaf)
+ : local_workspace(gaf.local_workspace), expr(gaf.expr), gis(0)
+ { if (gaf.gis) compile(); }
+
+ void ga_function::compile() const {
+ if (gis) delete gis;
+ gis = new ga_instruction_set;
+ local_workspace.clear_expressions();
+ local_workspace.add_function_expression(expr);
+ ga_compile_function(local_workspace, *gis, false);
+ }
+
+ ga_function &ga_function::operator =(const ga_function &gaf) {
+ if (gis) delete gis;
+ gis = 0;
+ local_workspace = gaf.local_workspace;
+ expr = gaf.expr;
+ if (gaf.gis) compile();
+ return *this;
+ }
+
+ ga_function::~ga_function() { if (gis) delete gis; gis = 0; }
+
+ const base_tensor &ga_function::eval() const {
+ GMM_ASSERT1(gis, "Uncompiled function");
+ gmm::clear(local_workspace.assembled_tensor().as_vector());
+ ga_function_exec(*gis);
+ return local_workspace.assembled_tensor();
+ }
+
+ void ga_function::derivative(const std::string &var) {
+ GMM_ASSERT1(gis, "Uncompiled function");
+ if (local_workspace.nb_trees()) {
+ ga_tree tree = *(local_workspace.tree_info(0).ptree);
+ ga_derivative(tree, local_workspace, *((const mesh *)(0)), var, "", 1);
+ if (tree.root) {
+ ga_semantic_analysis(expr, tree, local_workspace, 1, 1, false, true);
+ // To be improved to suppress test functions in the expression ...
+ // ga_replace_test_by_cte do not work in all operations like
+ // vector components x(1)
+ // ga_replace_test_by_cte(tree.root, false);
+ // ga_semantic_analysis(expr, tree, local_workspace, 1, 1,
+ // false, true);
+ }
+ expr = ga_tree_to_string(tree);
+ }
+ if (gis) delete gis;
+ gis = 0;
+ compile();
+ }
+
} /* end of namespace */
diff --git a/src/getfem_generic_assembly_interpolation.cc
b/src/getfem_generic_assembly_interpolation.cc
new file mode 100644
index 0000000..b6bde52
--- /dev/null
+++ b/src/getfem_generic_assembly_interpolation.cc
@@ -0,0 +1,914 @@
+/*===========================================================================
+
+ Copyright (C) 2013-2018 Yves Renard
+
+ This file is a part of GetFEM++
+
+ GetFEM++ is free software; you can redistribute it and/or modify it
+ under the terms of the GNU Lesser General Public License as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version along with the GCC Runtime Library
+ Exception either version 3.1 or (at your option) any later version.
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+ License and GCC Runtime Library Exception for more details.
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+===========================================================================*/
+
+#include "getfem/getfem_generic_assembly_tree.h"
+#include "getfem/getfem_generic_assembly_semantic.h"
+#include "getfem/getfem_generic_assembly_compile_and_exec.h"
+#include "getfem/getfem_generic_assembly_functions_and_operators.h"
+
+namespace getfem {
+
+ //=========================================================================
+ // Interpolation functions
+ //=========================================================================
+
+ // general Interpolation
+ void ga_interpolation(ga_workspace &workspace,
+ ga_interpolation_context &gic) {
+ ga_instruction_set gis;
+ ga_compile_interpolation(workspace, gis);
+ ga_interpolation_exec(gis, workspace, gic);
+ }
+
+ // Interpolation on a Lagrange fem on the same mesh
+ struct ga_interpolation_context_fem_same_mesh
+ : public ga_interpolation_context {
+ base_vector &result;
+ std::vector<int> dof_count;
+ const mesh_fem &mf;
+ bool initialized;
+ bool is_torus;
+ size_type s;
+
+ virtual bgeot::pstored_point_tab
+ ppoints_for_element(size_type cv, short_type f,
+ std::vector<size_type> &ind) const {
+ pfem pf = mf.fem_of_element(cv);
+ GMM_ASSERT1(pf->is_lagrange(),
+ "Only Lagrange fems can be used in interpolation");
+
+ if (f != short_type(-1)) {
+
+ for (size_type i = 0;
+ i < pf->node_convex(cv).structure()->nb_points_of_face(f); ++i)
+ ind.push_back
+ (pf->node_convex(cv).structure()->ind_points_of_face(f)[i]);
+ } else {
+ for (size_type i = 0; i < pf->node_convex(cv).nb_points(); ++i)
+ ind.push_back(i);
+ }
+
+ return pf->node_tab(cv);
+ }
+
+ virtual bool use_pgp(size_type) const { return true; }
+ virtual bool use_mim() const { return false; }
+
+ void init_(size_type si, size_type q, size_type qmult) {
+ s = si;
+ gmm::resize(result, qmult * mf.nb_basic_dof());
+ gmm::clear(result);
+ gmm::resize(dof_count, mf.nb_basic_dof()/q);
+ gmm::clear(dof_count);
+ initialized = true;
+ }
+
+ void store_result_for_torus(size_type cv, size_type i, base_tensor &t) {
+ size_type target_dim = mf.fem_of_element(cv)->target_dim();
+ GMM_ASSERT2(target_dim == 3, "Invalid torus fem.");
+ size_type qdim = 1;
+ size_type result_dim = 2;
+ if (!initialized) {init_(qdim, qdim, qdim);}
+ size_type idof = mf.ind_basic_dof_of_element(cv)[i];
+ result[idof] = t[idof%result_dim];
+ ++dof_count[idof];
+ }
+
+ virtual void store_result(size_type cv, size_type i, base_tensor &t) {
+ if (is_torus){
+ store_result_for_torus(cv, i, t);
+ return;
+ }
+ size_type si = t.size();
+ size_type q = mf.get_qdim();
+ size_type qmult = si / q;
+ GMM_ASSERT1( (si % q) == 0, "Incompatibility between the mesh_fem and "
+ "the size of the expression to be interpolated");
+ if (!initialized) { init_(si, q, qmult); }
+ GMM_ASSERT1(s == si, "Internal error");
+ size_type idof = mf.ind_basic_dof_of_element(cv)[i*q];
+ gmm::add(t.as_vector(),
+ gmm::sub_vector(result, gmm::sub_interval(qmult*idof, s)));
+ (dof_count[idof/q])++;
+ }
+
+ virtual void finalize() {
+ std::vector<size_type> data(3);
+ data[0] = initialized ? result.size() : 0;
+ data[1] = initialized ? dof_count.size() : 0;
+ data[2] = initialized ? s : 0;
+ MPI_MAX_VECTOR(data);
+ if (!initialized) {
+ if (data[0]) {
+ gmm::resize(result, data[0]);
+ gmm::resize(dof_count, data[1]);
+ gmm::clear(dof_count);
+ s = data[2];
+ }
+ gmm::clear(result);
+ }
+ if (initialized)
+ GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[2] &&
+ gmm::vect_size(dof_count) == data[1], "Incompatible sizes");
+ MPI_SUM_VECTOR(result);
+ MPI_SUM_VECTOR(dof_count);
+ for (size_type i = 0; i < dof_count.size(); ++i)
+ if (dof_count[i])
+ gmm::scale(gmm::sub_vector(result, gmm::sub_interval(s*i, s)),
+ scalar_type(1) / scalar_type(dof_count[i]));
+ }
+
+ virtual const mesh &linked_mesh() { return mf.linked_mesh(); }
+
+ ga_interpolation_context_fem_same_mesh(const mesh_fem &mf_, base_vector &r)
+ : result(r), mf(mf_), initialized(false), is_torus(false) {
+ GMM_ASSERT1(!(mf.is_reduced()),
+ "Interpolation on reduced fem is not allowed");
+ if (dynamic_cast<const torus_mesh_fem*>(&mf)){
+ auto first_cv = mf.first_convex_of_basic_dof(0);
+ auto target_dim = mf.fem_of_element(first_cv)->target_dim();
+ if (target_dim == 3) is_torus = true;
+ }
+ }
+ };
+
+ void ga_interpolation_Lagrange_fem
+ (ga_workspace &workspace, const mesh_fem &mf, base_vector &result) {
+ ga_interpolation_context_fem_same_mesh gic(mf, result);
+ ga_interpolation(workspace, gic);
+ }
+
+ void ga_interpolation_Lagrange_fem
+ (const getfem::model &md, const std::string &expr, const mesh_fem &mf,
+ base_vector &result, const mesh_region &rg) {
+ ga_workspace workspace(md);
+ workspace.add_interpolation_expression(expr, mf.linked_mesh(), rg);
+ ga_interpolation_Lagrange_fem(workspace, mf, result);
+ }
+
+ // Interpolation on a cloud of points
+ struct ga_interpolation_context_mti
+ : public ga_interpolation_context {
+ base_vector &result;
+ const mesh_trans_inv &mti;
+ bool initialized;
+ size_type s, nbdof;
+
+
+ virtual bgeot::pstored_point_tab
+ ppoints_for_element(size_type cv, short_type,
+ std::vector<size_type> &ind) const {
+ std::vector<size_type> itab;
+ mti.points_on_convex(cv, itab);
+ std::vector<base_node> pt_tab(itab.size());
+ for (size_type i = 0; i < itab.size(); ++i) {
+ pt_tab[i] = mti.reference_coords()[itab[i]];
+ ind.push_back(i);
+ }
+ return store_point_tab(pt_tab);
+ }
+
+ virtual bool use_pgp(size_type) const { return false; }
+ virtual bool use_mim() const { return false; }
+
+ virtual void store_result(size_type cv, size_type i, base_tensor &t) {
+ size_type si = t.size();
+ if (!initialized) {
+ s = si;
+ gmm::resize(result, s * nbdof);
+ gmm::clear(result);
+ initialized = true;
+ }
+ GMM_ASSERT1(s == si, "Internal error");
+ size_type ipt = mti.point_on_convex(cv, i);
+ size_type dof_t = mti.id_of_point(ipt);
+ gmm::add(t.as_vector(),
+ gmm::sub_vector(result, gmm::sub_interval(s*dof_t, s)));
+ }
+
+ virtual void finalize() {
+ std::vector<size_type> data(2);
+ data[0] = initialized ? result.size() : 0;
+ data[1] = initialized ? s : 0;
+ MPI_MAX_VECTOR(data);
+ if (!initialized) {
+ if (data[0]) {
+ gmm::resize(result, data[0]);
+ s = data[1];
+ }
+ gmm::clear(result);
+ }
+ if (initialized)
+ GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[1],
+ "Incompatible sizes");
+ MPI_SUM_VECTOR(result);
+ }
+
+ virtual const mesh &linked_mesh() { return mti.linked_mesh(); }
+
+ ga_interpolation_context_mti(const mesh_trans_inv &mti_, base_vector &r,
+ size_type nbdof_ = size_type(-1))
+ : result(r), mti(mti_), initialized(false), nbdof(nbdof_) {
+ if (nbdof == size_type(-1)) nbdof = mti.nb_points();
+ }
+ };
+
+ // Distribute to be parallelized
+ void ga_interpolation_mti
+ (const getfem::model &md, const std::string &expr, mesh_trans_inv &mti,
+ base_vector &result, const mesh_region &rg, int extrapolation,
+ const mesh_region &rg_source, size_type nbdof) {
+
+ ga_workspace workspace(md);
+ workspace.add_interpolation_expression(expr, mti.linked_mesh(), rg);
+
+ mti.distribute(extrapolation, rg_source);
+ ga_interpolation_context_mti gic(mti, result, nbdof);
+ ga_interpolation(workspace, gic);
+ }
+
+
+ // Interpolation on a im_data
+ struct ga_interpolation_context_im_data
+ : public ga_interpolation_context {
+ base_vector &result;
+ const im_data &imd;
+ bool initialized;
+ size_type s;
+
+ virtual bgeot::pstored_point_tab
+ ppoints_for_element(size_type cv, short_type f,
+ std::vector<size_type> &ind) const {
+ pintegration_method pim =imd.linked_mesh_im().int_method_of_element(cv);
+ if (pim->type() == IM_NONE) return bgeot::pstored_point_tab();
+ GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot "
+ "be used in high level generic assembly");
+ size_type i_start(0), i_end(0);
+ if (f == short_type(-1))
+ i_end = pim->approx_method()->nb_points_on_convex();
+ else {
+ i_start = pim->approx_method()->ind_first_point_on_face(f);
+ i_end = i_start + pim->approx_method()->nb_points_on_face(f);
+ }
+ for (size_type i = i_start; i < i_end; ++i) ind.push_back(i);
+ return pim->approx_method()->pintegration_points();
+ }
+
+ virtual bool use_pgp(size_type cv) const {
+ pintegration_method pim =imd.linked_mesh_im().int_method_of_element(cv);
+ if (pim->type() == IM_NONE) return false;
+ GMM_ASSERT1(pim->type() == IM_APPROX, "Sorry, exact methods cannot "
+ "be used in high level generic assembly");
+ return !(pim->approx_method()->is_built_on_the_fly());
+ }
+ virtual bool use_mim() const { return true; }
+
+ virtual void store_result(size_type cv, size_type i, base_tensor &t) {
+ size_type si = t.size();
+ if (!initialized) {
+ s = si;
+ GMM_ASSERT1(imd.tensor_size() == t.sizes() ||
+ (imd.tensor_size().size() == size_type(1) &&
+ imd.tensor_size()[0] == size_type(1) &&
+ si == size_type(1)),
+ "Im_data tensor size " << imd.tensor_size() <<
+ " does not match the size of the interpolated "
+ "expression " << t.sizes() << ".");
+ gmm::resize(result, s * imd.nb_filtered_index());
+ gmm::clear(result);
+ initialized = true;
+ }
+ GMM_ASSERT1(s == si, "Internal error");
+ size_type ipt = imd.filtered_index_of_point(cv, i);
+ GMM_ASSERT1(ipt != size_type(-1),
+ "Im data with no data on the current integration point.");
+ gmm::add(t.as_vector(),
+ gmm::sub_vector(result, gmm::sub_interval(s*ipt, s)));
+ }
+
+ virtual void finalize() {
+ std::vector<size_type> data(2);
+ data[0] = initialized ? result.size() : 0;
+ data[1] = initialized ? s : 0;
+ MPI_MAX_VECTOR(data);
+ if (initialized) {
+ GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[1],
+ "Incompatible sizes");
+ } else {
+ if (data[0]) {
+ gmm::resize(result, data[0]);
+ s = data[1];
+ }
+ gmm::clear(result);
+ }
+ MPI_SUM_VECTOR(result);
+ }
+
+ virtual const mesh &linked_mesh() { return imd.linked_mesh(); }
+
+ ga_interpolation_context_im_data(const im_data &imd_, base_vector &r)
+ : result(r), imd(imd_), initialized(false) { }
+ };
+
+ void ga_interpolation_im_data
+ (ga_workspace &workspace, const im_data &imd, base_vector &result) {
+ ga_interpolation_context_im_data gic(imd, result);
+ ga_interpolation(workspace, gic);
+ }
+
+ void ga_interpolation_im_data
+ (const getfem::model &md, const std::string &expr, const im_data &imd,
+ base_vector &result, const mesh_region &rg) {
+ ga_workspace workspace(md);
+ workspace.add_interpolation_expression
+ (expr, imd.linked_mesh_im(), rg);
+
+ ga_interpolation_im_data(workspace, imd, result);
+ }
+
+
+ // Interpolation on a stored_mesh_slice
+ struct ga_interpolation_context_mesh_slice
+ : public ga_interpolation_context {
+ base_vector &result;
+ const stored_mesh_slice &sl;
+ bool initialized;
+ size_type s;
+ std::vector<size_type> first_node;
+
+ virtual bgeot::pstored_point_tab
+ ppoints_for_element(size_type cv, short_type f,
+ std::vector<size_type> &ind) const {
+ GMM_ASSERT1(f == short_type(-1), "No support for interpolation on faces"
+ " for a stored_mesh_slice yet.");
+ size_type ic = sl.convex_pos(cv);
+ const mesh_slicer::cs_nodes_ct &nodes = sl.nodes(ic);
+ std::vector<base_node> pt_tab(nodes.size());
+ for (size_type i=0; i < nodes.size(); ++i) {
+ pt_tab[i] = nodes[i].pt_ref;
+ ind.push_back(i);
+ }
+ return store_point_tab(pt_tab);
+ }
+
+ virtual bool use_pgp(size_type /* cv */) const { return false; } // why
not?
+ virtual bool use_mim() const { return false; }
+
+ virtual void store_result(size_type cv, size_type i, base_tensor &t) {
+ size_type si = t.size();
+ if (!initialized) {
+ s = si;
+ gmm::resize(result, s * sl.nb_points());
+ gmm::clear(result);
+ initialized = true;
+ first_node.resize(sl.nb_convex());
+ for (size_type ic=0; ic < sl.nb_convex()-1; ++ic)
+ first_node[ic+1] = first_node[ic] + sl.nodes(ic).size();
+ }
+ GMM_ASSERT1(s == si && result.size() == s * sl.nb_points(), "Internal
error");
+ size_type ic = sl.convex_pos(cv);
+ size_type ipt = first_node[ic] + i;
+ gmm::add(t.as_vector(),
+ gmm::sub_vector(result, gmm::sub_interval(s*ipt, s)));
+ }
+
+ virtual void finalize() {
+ std::vector<size_type> data(2);
+ data[0] = initialized ? result.size() : 0;
+ data[1] = initialized ? s : 0;
+ MPI_MAX_VECTOR(data);
+ if (initialized) {
+ GMM_ASSERT1(gmm::vect_size(result) == data[0] && s == data[1],
+ "Incompatible sizes");
+ } else {
+ if (data[0]) {
+ gmm::resize(result, data[0]);
+ s = data[1];
+ }
+ gmm::clear(result);
+ }
+ MPI_SUM_VECTOR(result);
+ }
+
+ virtual const mesh &linked_mesh() { return sl.linked_mesh(); }
+
+ ga_interpolation_context_mesh_slice(const stored_mesh_slice &sl_,
base_vector &r)
+ : result(r), sl(sl_), initialized(false) { }
+ };
+
+ void ga_interpolation_mesh_slice
+ (ga_workspace &workspace, const stored_mesh_slice &sl, base_vector &result) {
+ ga_interpolation_context_mesh_slice gic(sl, result);
+ ga_interpolation(workspace, gic);
+ }
+
+ void ga_interpolation_mesh_slice
+ (const getfem::model &md, const std::string &expr, const stored_mesh_slice
&sl,
+ base_vector &result, const mesh_region &rg) {
+ ga_workspace workspace(md);
+ workspace.add_interpolation_expression(expr, sl.linked_mesh(), rg);
+ ga_interpolation_mesh_slice(workspace, sl, result);
+ }
+
+
+ //=========================================================================
+ // Local projection functions
+ //=========================================================================
+
+ void ga_local_projection(const getfem::model &md, const mesh_im &mim,
+ const std::string &expr, const mesh_fem &mf,
+ base_vector &result, const mesh_region ®ion) {
+
+ // The case where the expression is a vector one and mf a scalar fem is
+ // not taken into account for the moment.
+
+ // Could be improved by not performing the assembly of the global mass
matrix
+ // working locally. This means a specific assembly.
+ model_real_sparse_matrix M(mf.nb_dof(), mf.nb_dof());
+ asm_mass_matrix(M, mim, mf, region);
+
+ ga_workspace workspace(md);
+ size_type nbdof = md.nb_dof();
+ gmm::sub_interval I(nbdof, mf.nb_dof());
+ workspace.add_fem_variable("c__dummy_var_95_", mf, I, base_vector(nbdof));
+ if (mf.get_qdims().size() > 1)
+
workspace.add_expression("("+expr+"):Test_c__dummy_var_95_",mim,region,2);
+ else
+
workspace.add_expression("("+expr+").Test_c__dummy_var_95_",mim,region,2);
+ base_vector residual(nbdof+mf.nb_dof());
+ workspace.set_assembled_vector(residual);
+ workspace.assembly(1);
+ getfem::base_vector F(mf.nb_dof());
+ gmm::resize(result, mf.nb_dof());
+ gmm::copy(gmm::sub_vector(residual, I), F);
+
+ getfem::base_matrix loc_M;
+ getfem::base_vector loc_U;
+ for (mr_visitor v(region, mf.linked_mesh(), true); !v.finished(); ++v) {
+ if (mf.convex_index().is_in(v.cv())) {
+ size_type nd = mf.nb_basic_dof_of_element(v.cv());
+ loc_M.base_resize(nd, nd); gmm::resize(loc_U, nd);
+ gmm::sub_index J(mf.ind_basic_dof_of_element(v.cv()));
+ gmm::copy(gmm::sub_matrix(M, J, J), loc_M);
+ gmm::lu_solve(loc_M, loc_U, gmm::sub_vector(F, J));
+ gmm::copy(loc_U, gmm::sub_vector(result, J));
+ }
+ }
+ MPI_SUM_VECTOR(result);
+ }
+
+ //=========================================================================
+ // Interpolate transformation with an expression
+ //=========================================================================
+
+ class interpolate_transformation_expression
+ : public virtual_interpolate_transformation, public context_dependencies {
+
+ struct workspace_gis_pair : public std::pair<ga_workspace,
ga_instruction_set> {
+ inline ga_workspace &workspace() { return this->first; }
+ inline ga_instruction_set &gis() { return this->second; }
+ };
+
+ const mesh &source_mesh;
+ const mesh &target_mesh;
+ std::string expr;
+ mutable bgeot::rtree element_boxes;
+ mutable bool recompute_elt_boxes;
+ mutable ga_workspace local_workspace;
+ mutable ga_instruction_set local_gis;
+ mutable bgeot::geotrans_inv_convex gic;
+ mutable base_node P;
+ mutable std::set<var_trans_pair> used_vars;
+ mutable std::set<var_trans_pair> used_data;
+ mutable std::map<var_trans_pair,
+ workspace_gis_pair> compiled_derivatives;
+ mutable bool extract_variable_done;
+ mutable bool extract_data_done;
+
+ public:
+ void update_from_context() const {
+ recompute_elt_boxes = true;
+ }
+
+ void extract_variables(const ga_workspace &workspace,
+ std::set<var_trans_pair> &vars,
+ bool ignore_data, const mesh &/* m */,
+ const std::string &/* interpolate_name */) const {
+ if ((ignore_data && !extract_variable_done) ||
+ (!ignore_data && !extract_data_done)) {
+ if (ignore_data)
+ used_vars.clear();
+ else
+ used_data.clear();
+ ga_workspace aux_workspace;
+ aux_workspace = ga_workspace(true, workspace);
+ aux_workspace.clear_expressions();
+ aux_workspace.add_interpolation_expression(expr, source_mesh);
+ for (size_type i = 0; i < aux_workspace.nb_trees(); ++i)
+ ga_extract_variables(aux_workspace.tree_info(i).ptree->root,
+ aux_workspace, source_mesh,
+ ignore_data ? used_vars : used_data,
+ ignore_data);
+ if (ignore_data)
+ extract_variable_done = true;
+ else
+ extract_data_done = true;
+ }
+ if (ignore_data)
+ vars.insert(used_vars.begin(), used_vars.end());
+ else
+ vars.insert(used_data.begin(), used_data.end());
+ }
+
+ void init(const ga_workspace &workspace) const {
+ size_type N = target_mesh.dim();
+
+ // Expression compilation
+ local_workspace = ga_workspace(true, workspace);
+ local_workspace.clear_expressions();
+
+ local_workspace.add_interpolation_expression(expr, source_mesh);
+ local_gis = ga_instruction_set();
+ ga_compile_interpolation(local_workspace, local_gis);
+
+ // In fact, transformations are not allowed ... for future compatibility
+ for (const std::string &transname : local_gis.transformations)
+ local_workspace.interpolate_transformation(transname)
+ ->init(local_workspace);
+
+ if (!extract_variable_done) {
+ std::set<var_trans_pair> vars;
+ extract_variables(workspace, vars, true, source_mesh, "");
+ }
+
+ for (const var_trans_pair &var : used_vars) {
+ workspace_gis_pair &pwi = compiled_derivatives[var];
+ pwi.workspace() = local_workspace;
+ pwi.gis() = ga_instruction_set();
+ if (pwi.workspace().nb_trees()) {
+ ga_tree tree = *(pwi.workspace().tree_info(0).ptree);
+ ga_derivative(tree, pwi.workspace(), source_mesh,
+ var.varname, var.transname, 1);
+ if (tree.root)
+ ga_semantic_analysis(expr, tree, local_workspace, 1, 1,
+ false, true);
+ ga_compile_interpolation(pwi.workspace(), pwi.gis());
+ }
+ }
+
+ // Element_boxes update (if necessary)
+ if (recompute_elt_boxes) {
+
+ element_boxes.clear();
+ base_node bmin(N), bmax(N);
+ for (dal::bv_visitor cv(target_mesh.convex_index());
+ !cv.finished(); ++cv) {
+
+ bgeot::pgeometric_trans pgt = target_mesh.trans_of_convex(cv);
+
+ size_type nbd_t = pgt->nb_points();
+ if (nbd_t) {
+ gmm::copy(target_mesh.points_of_convex(cv)[0], bmin);
+ gmm::copy(bmin, bmax);
+ } else {
+ gmm::clear(bmin);
+ gmm::clear(bmax);
+ }
+ for (short_type ip = 1; ip < nbd_t; ++ip) {
+ // size_type ind = target_mesh.ind_points_of_convex(cv)[ip];
+ const base_node &pt = target_mesh.points_of_convex(cv)[ip];
+
+ for (size_type k = 0; k < N; ++k) {
+ bmin[k] = std::min(bmin[k], pt[k]);
+ bmax[k] = std::max(bmax[k], pt[k]);
+ }
+ }
+
+ scalar_type h = bmax[0] - bmin[0];
+ for (size_type k = 1; k < N; ++k) h = std::max(h, bmax[k]-bmin[k]);
+ if (pgt->is_linear()) h *= 1E-4;
+ for (auto &&val : bmin) val -= h*0.2;
+ for (auto &&val : bmax) val += h*0.2;
+
+ element_boxes.add_box(bmin, bmax, cv);
+ }
+ recompute_elt_boxes = false;
+ }
+ }
+
+ void finalize() const {
+ for (const std::string &transname : local_gis.transformations)
+ local_workspace.interpolate_transformation(transname)->finalize();
+ local_gis = ga_instruction_set();
+ }
+
+
+ int transform(const ga_workspace &/*workspace*/, const mesh &m,
+ fem_interpolation_context &ctx_x,
+ const base_small_vector &Normal,
+ const mesh **m_t,
+ size_type &cv, short_type &face_num,
+ base_node &P_ref,
+ base_small_vector &/*N_y*/,
+ std::map<var_trans_pair, base_tensor> &derivatives,
+ bool compute_derivatives) const {
+ int ret_type = 0;
+
+ ga_interpolation_single_point_exec(local_gis, local_workspace, ctx_x,
+ Normal, m);
+
+ GMM_ASSERT1(local_workspace.assembled_tensor().size() == m.dim(),
+ "Wrong dimension of the tranformation expression");
+ P.resize(m.dim());
+ gmm::copy(local_workspace.assembled_tensor().as_vector(), P);
+
+ bgeot::rtree::pbox_set bset;
+ element_boxes.find_boxes_at_point(P, bset);
+ *m_t = &target_mesh;
+
+ while (bset.size()) {
+ bgeot::rtree::pbox_set::iterator it = bset.begin(), itmax = it;
+
+ if (bset.size() > 1) {
+ // Searching the box for which the point is the most in the interior
+ scalar_type rate_max = scalar_type(-1);
+ for (; it != bset.end(); ++it) {
+
+ scalar_type rate_box = scalar_type(1);
+ for (size_type i = 0; i < m.dim(); ++i) {
+ scalar_type h = (*it)->max[i] - (*it)->min[i];
+ if (h > scalar_type(0)) {
+ scalar_type rate
+ = std::min((*it)->max[i] - P[i], P[i] - (*it)->min[i]) / h;
+ rate_box = std::min(rate, rate_box);
+ }
+ }
+ if (rate_box > rate_max) {
+ itmax = it;
+ rate_max = rate_box;
+ }
+ }
+ }
+
+ cv = (*itmax)->id;
+ gic.init(target_mesh.points_of_convex(cv),
+ target_mesh.trans_of_convex(cv));
+
+ bool converged = true;
+ bool is_in = gic.invert(P, P_ref, converged, 1E-4);
+ // cout << "cv = " << cv << " P = " << P << " P_ref = " << P_ref <<
endl;
+ // cout << " is_in = " << int(is_in) << endl;
+ // for (size_type iii = 0;
+ // iii < target_mesh.points_of_convex(cv).size(); ++iii)
+ // cout << target_mesh.points_of_convex(cv)[iii] << endl;
+
+ if (is_in && converged) {
+ face_num = short_type(-1); // Should detect potential faces ?
+ ret_type = 1;
+ break;
+ }
+
+ if (bset.size() == 1) break;
+ bset.erase(itmax);
+ }
+
+ // Note on derivatives of the transformation : for efficiency and
+ // simplicity reasons, the derivative should be computed with
+ // the value of corresponding test functions. This means that
+ // for a transformation F(u) the computed derivative is F'(u).Test_u
+ // including the Test_u.
+ if (compute_derivatives) { // To be tested both with the computation
+ // of derivative. Could be optimized ?
+ for (auto &&d : derivatives) {
+ workspace_gis_pair &pwi = compiled_derivatives[d.first];
+
+ gmm::clear(pwi.workspace().assembled_tensor().as_vector());
+ ga_function_exec(pwi.gis());
+ d.second = pwi.workspace().assembled_tensor();
+ }
+ }
+ return ret_type;
+ }
+
+ interpolate_transformation_expression(const mesh &sm, const mesh &tm,
+ const std::string &expr_)
+ : source_mesh(sm), target_mesh(tm), expr(expr_),
+ recompute_elt_boxes(true), extract_variable_done(false),
+ extract_data_done(false)
+ { this->add_dependency(tm); }
+
+ };
+
+
+ void add_interpolate_transformation_from_expression
+ (model &md, const std::string &name, const mesh &sm, const mesh &tm,
+ const std::string &expr) {
+ pinterpolate_transformation
+ p = std::make_shared<interpolate_transformation_expression>(sm, tm,
expr);
+ md.add_interpolate_transformation(name, p);
+ }
+
+ void add_interpolate_transformation_from_expression
+ (ga_workspace &workspace, const std::string &name, const mesh &sm,
+ const mesh &tm, const std::string &expr) {
+ pinterpolate_transformation
+ p = std::make_shared<interpolate_transformation_expression>(sm, tm,
expr);
+ workspace.add_interpolate_transformation(name, p);
+ }
+
+ //=========================================================================
+ // Interpolate transformation on neighbour element (for internal faces)
+ //=========================================================================
+
+ class interpolate_transformation_neighbour
+ : public virtual_interpolate_transformation, public context_dependencies {
+
+ public:
+ void update_from_context() const {}
+ void extract_variables(const ga_workspace &/* workspace */,
+ std::set<var_trans_pair> &/* vars */,
+ bool /* ignore_data */, const mesh &/* m */,
+ const std::string &/* interpolate_name */) const {}
+ void init(const ga_workspace &/* workspace */) const {}
+ void finalize() const {}
+
+ int transform(const ga_workspace &/*workspace*/, const mesh &m_x,
+ fem_interpolation_context &ctx_x,
+ const base_small_vector &/*Normal*/, const mesh **m_t,
+ size_type &cv, short_type &face_num,
+ base_node &P_ref,
+ base_small_vector &/*N_y*/,
+ std::map<var_trans_pair, base_tensor> &/*derivatives*/,
+ bool compute_derivatives) const {
+
+ int ret_type = 0;
+ *m_t = &m_x;
+ size_type cv_x = ctx_x.convex_num();
+ short_type face_x = ctx_x.face_num();
+ GMM_ASSERT1(face_x != short_type(-1), "Neighbour transformation can "
+ "only be applied to internal faces");
+
+ auto adj_face = m_x.adjacent_face(cv_x, face_x);
+
+ if (adj_face.cv != size_type(-1)) {
+ bgeot::geotrans_inv_convex gic;
+ gic.init(m_x.points_of_convex(adj_face.cv),
+ m_x.trans_of_convex(adj_face.cv));
+ bool converged = true;
+ bool is_in = gic.invert(ctx_x.xreal(), P_ref, converged, 1E-4);
+ GMM_ASSERT1(is_in && converged, "Geometric transformation inversion "
+ "has failed in neighbour transformation");
+ face_num = adj_face.f;
+ cv = adj_face.cv;
+ ret_type = 1;
+ }
+ GMM_ASSERT1(!compute_derivatives,
+ "No derivative for this transformation");
+ return ret_type;
+ }
+
+ interpolate_transformation_neighbour() { }
+
+ };
+
+
+ pinterpolate_transformation interpolate_transformation_neighbour_instance()
+ {
+ pinterpolate_transformation
+ p = std::make_shared<interpolate_transformation_neighbour>();
+ return p;
+ }
+
+ //=========================================================================
+ // Interpolate transformation on neighbour element (for extrapolation)
+ //=========================================================================
+
+ class interpolate_transformation_element_extrapolation
+ : public virtual_interpolate_transformation, public context_dependencies {
+
+ const mesh &sm;
+ std::map<size_type, size_type> elt_corr;
+
+ public:
+ void update_from_context() const {}
+ void extract_variables(const ga_workspace &/* workspace */,
+ std::set<var_trans_pair> &/* vars */,
+ bool /* ignore_data */, const mesh &/* m */,
+ const std::string &/* interpolate_name */) const {}
+ void init(const ga_workspace &/* workspace */) const {}
+ void finalize() const {}
+
+ int transform(const ga_workspace &/*workspace*/, const mesh &m_x,
+ fem_interpolation_context &ctx_x,
+ const base_small_vector &/*Normal*/, const mesh **m_t,
+ size_type &cv, short_type &face_num,
+ base_node &P_ref,
+ base_small_vector &/*N_y*/,
+ std::map<var_trans_pair, base_tensor> &/*derivatives*/,
+ bool compute_derivatives) const {
+ int ret_type = 0;
+ *m_t = &m_x;
+ GMM_ASSERT1(&sm == &m_x, "Bad mesh");
+ size_type cv_x = ctx_x.convex_num(), cv_y = cv_x;
+ auto it = elt_corr.find(cv_x);
+ if (it != elt_corr.end()) cv_y = it->second;
+
+ if (cv_x != cv_y) {
+ bgeot::geotrans_inv_convex gic;
+ gic.init(m_x.points_of_convex(cv_y),
+ m_x.trans_of_convex(cv_y));
+ bool converged = true;
+ gic.invert(ctx_x.xreal(), P_ref, converged, 1E-4);
+ GMM_ASSERT1(converged, "Geometric transformation inversion "
+ "has failed in element extrapolation transformation");
+ face_num = short_type(-1);
+ cv = cv_y;
+ ret_type = 1;
+ } else {
+ cv = cv_x;
+ face_num = short_type(-1);
+ P_ref = ctx_x.xref();
+ ret_type = 1;
+ }
+ GMM_ASSERT1(!compute_derivatives,
+ "No derivative for this transformation");
+ return ret_type;
+ }
+
+ void set_correspondance(const std::map<size_type, size_type> &ec) {
+ elt_corr = ec;
+ }
+
+ interpolate_transformation_element_extrapolation
+ (const mesh &sm_, const std::map<size_type, size_type> &ec)
+ : sm(sm_), elt_corr(ec) { }
+ };
+
+
+ void add_element_extrapolation_transformation
+ (model &md, const std::string &name, const mesh &sm,
+ std::map<size_type, size_type> &elt_corr) {
+ pinterpolate_transformation
+ p = std::make_shared<interpolate_transformation_element_extrapolation>
+ (sm, elt_corr);
+ md.add_interpolate_transformation(name, p);
+ }
+
+ void add_element_extrapolation_transformation
+ (ga_workspace &workspace, const std::string &name, const mesh &sm,
+ std::map<size_type, size_type> &elt_corr) {
+ pinterpolate_transformation
+ p = std::make_shared<interpolate_transformation_element_extrapolation>
+ (sm, elt_corr);
+ workspace.add_interpolate_transformation(name, p);
+ }
+
+ void set_element_extrapolation_correspondance
+ (ga_workspace &workspace, const std::string &name,
+ std::map<size_type, size_type> &elt_corr) {
+ GMM_ASSERT1(workspace.interpolate_transformation_exists(name),
+ "Unknown transformation");
+ auto pit=workspace.interpolate_transformation(name).get();
+ auto cpext
+ = dynamic_cast<const interpolate_transformation_element_extrapolation *>
+ (pit);
+ GMM_ASSERT1(cpext,
+ "The transformation is not of element extrapolation type");
+ const_cast<interpolate_transformation_element_extrapolation *>(cpext)
+ ->set_correspondance(elt_corr);
+ }
+
+ void set_element_extrapolation_correspondance
+ (model &md, const std::string &name,
+ std::map<size_type, size_type> &elt_corr) {
+ GMM_ASSERT1(md.interpolate_transformation_exists(name),
+ "Unknown transformation");
+ auto pit=md.interpolate_transformation(name).get();
+ auto cpext
+ = dynamic_cast<const interpolate_transformation_element_extrapolation *>
+ (pit);
+ GMM_ASSERT1(cpext,
+ "The transformation is not of element extrapolation type");
+ const_cast<interpolate_transformation_element_extrapolation *>(cpext)
+ ->set_correspondance(elt_corr);
+ }
+
+} /* end of namespace */
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index 006d573..787337c 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -24,10 +24,14 @@
#include <getfem/getfem_generic_assembly_functions_and_operators.h>
#include <getfem/getfem_generic_assembly_semantic.h>
+#include <getfem/getfem_generic_assembly_compile_and_exec.h>
namespace getfem {
-
+ extern bool predef_operators_nonlinear_elasticity_initialized;
+ extern bool predef_operators_plasticity_initialized;
+ extern bool predef_operators_contact_initialized;
+
bool ga_extract_variables(const pga_tree_node pnode,
const ga_workspace &workspace,
const mesh &m,
@@ -209,7 +213,9 @@ namespace getfem {
= dal::singleton<ga_predef_function_tab>::instance(0);
const ga_predef_operator_tab &PREDEF_OPERATORS
= dal::singleton<ga_predef_operator_tab>::instance(0);
-
+ const ga_spec_function_tab &SPEC_FUNCTIONS
+ = dal::singleton<ga_spec_function_tab>::instance(0);
+
switch (pnode->node_type) {
case GA_NODE_PREDEF_FUNC: case GA_NODE_OPERATOR: case GA_NODE_SPEC_FUNC :
case GA_NODE_CONSTANT: case GA_NODE_X: case GA_NODE_ELT_SIZE:
@@ -1930,10 +1936,6 @@ namespace getfem {
// cout << " final hash code = " << pnode->hash_value << endl;
}
- extern bool predef_operators_nonlinear_elasticity_initialized;
- extern bool predef_operators_plasticity_initialized;
- extern bool predef_operators_contact_initialized;
- extern bool predef_functions_initialized;
void ga_semantic_analysis(const std::string &expr, ga_tree &tree,
const ga_workspace &workspace,
@@ -1941,8 +1943,7 @@ namespace getfem {
size_type ref_elt_dim,
bool eval_fixed_size,
bool ignore_X, int option) {
- GMM_ASSERT1(predef_functions_initialized &&
- predef_operators_nonlinear_elasticity_initialized &&
+ GMM_ASSERT1(predef_operators_nonlinear_elasticity_initialized &&
predef_operators_plasticity_initialized &&
predef_operators_contact_initialized, "Internal error");
if (!(tree.root)) return;
diff --git a/src/getfem_generic_assembly_tree.cc
b/src/getfem_generic_assembly_tree.cc
index 7f9ce6c..99e6307 100644
--- a/src/getfem_generic_assembly_tree.cc
+++ b/src/getfem_generic_assembly_tree.cc
@@ -21,6 +21,7 @@
#include "getfem/getfem_generic_assembly_tree.h"
#include "getfem/getfem_generic_assembly_functions_and_operators.h"
+#include "getfem/getfem_generic_assembly_compile_and_exec.h"
namespace getfem {
@@ -1115,14 +1116,17 @@ namespace getfem {
if (name.compare(0, 11, "Derivative_") == 0)
return 2;
- ga_predef_function_tab &PREDEF_FUNCTIONS
+ const ga_predef_function_tab &PREDEF_FUNCTIONS
= dal::singleton<ga_predef_function_tab>::instance(0);
- ga_predef_operator_tab &PREDEF_OPERATORS
+ const ga_predef_operator_tab &PREDEF_OPERATORS
= dal::singleton<ga_predef_operator_tab>::instance(0);
+ const ga_spec_function_tab &SPEC_FUNCTIONS
+ = dal::singleton<ga_spec_function_tab>::instance(0);
+
ga_predef_function_tab::const_iterator it=PREDEF_FUNCTIONS.find(name);
if (it != PREDEF_FUNCTIONS.end())
return 1;
-
+
if (SPEC_FUNCTIONS.find(name) != SPEC_FUNCTIONS.end())
return 1;
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
new file mode 100644
index 0000000..ccccb9d
--- /dev/null
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -0,0 +1,1005 @@
+/*===========================================================================
+
+ Copyright (C) 2013-2018 Yves Renard
+
+ This file is a part of GetFEM++
+
+ GetFEM++ is free software; you can redistribute it and/or modify it
+ under the terms of the GNU Lesser General Public License as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version along with the GCC Runtime Library
+ Exception either version 3.1 or (at your option) any later version.
+ This program is distributed in the hope that it will be useful, but
+ WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+ or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+ License and GCC Runtime Library Exception for more details.
+ You should have received a copy of the GNU Lesser General Public License
+ along with this program; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+
+===========================================================================*/
+
+#include "getfem/getfem_generic_assembly_tree.h"
+#include "getfem/getfem_generic_assembly_semantic.h"
+#include "getfem/getfem_generic_assembly_compile_and_exec.h"
+#include "getfem/getfem_generic_assembly_functions_and_operators.h"
+
+namespace getfem {
+
+ //=========================================================================
+ // Structure dealing with user defined environment : constant, variables,
+ // functions, operators.
+ //=========================================================================
+
+ void ga_workspace::init() {
+ // allocate own storage for K an V to be used unless/until external
+ // storage is provided with set_assembled_matrix/vector
+ K = std::make_shared<model_real_sparse_matrix>(2,2);
+ V = std::make_shared<base_vector>(2);
+ // add default transformations
+ add_interpolate_transformation
+ ("neighbour_elt", interpolate_transformation_neighbour_instance());
+ }
+
+ // variables and variable groups
+ 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) {
+ variables[name] = var_description(true, true, &mf, I, &VV, 0, 1);
+ }
+
+ void ga_workspace::add_fixed_size_variable
+ (const std::string &name,
+ const gmm::sub_interval &I, const model_real_plain_vector &VV) {
+ variables[name] = var_description(true, false, 0, I, &VV, 0,
+ dim_type(gmm::vect_size(VV)));
+ }
+
+ void ga_workspace::add_fem_constant
+ (const std::string &name, const mesh_fem &mf,
+ const model_real_plain_vector &VV) {
+ GMM_ASSERT1(mf.nb_dof(), "The provided mesh_fem of variable" << name
+ << "has zero degrees of freedom.");
+ size_type Q = gmm::vect_size(VV)/mf.nb_dof();
+ if (Q == 0) Q = size_type(1);
+ variables[name] = var_description(false, true, &mf,
+ gmm::sub_interval(), &VV, 0, Q);
+ }
+
+ void ga_workspace::add_fixed_size_constant
+ (const std::string &name, const model_real_plain_vector &VV) {
+ variables[name] = var_description(false, false, 0,
+ gmm::sub_interval(), &VV, 0,
+ gmm::vect_size(VV));
+ }
+
+ void ga_workspace::add_im_data(const std::string &name, const im_data &imd,
+ const model_real_plain_vector &VV) {
+ variables[name] = var_description
+ (false, false, 0, gmm::sub_interval(), &VV, &imd,
+ gmm::vect_size(VV)/(imd.nb_filtered_index() * imd.nb_tensor_elem()));
+ }
+
+
+ bool ga_workspace::variable_exists(const std::string &name) const {
+ return (md && md->variable_exists(name)) ||
+ (parent_workspace && parent_workspace->variable_exists(name)) ||
+ (variables.find(name) != variables.end());
+ }
+
+ bool ga_workspace::variable_group_exists(const std::string &name) const {
+ return (variable_groups.find(name) != variable_groups.end()) ||
+ (md && md->variable_group_exists(name)) ||
+ (parent_workspace && parent_workspace->variable_group_exists(name));
+ }
+
+ const std::vector<std::string>&
+ ga_workspace::variable_group(const std::string &group_name) const {
+ std::map<std::string, std::vector<std::string> >::const_iterator
+ it = variable_groups.find(group_name);
+ if (it != variable_groups.end())
+ return (variable_groups.find(group_name))->second;
+ if (md && md->variable_group_exists(group_name))
+ return md->variable_group(group_name);
+ if (parent_workspace &&
+ parent_workspace->variable_group_exists(group_name))
+ return parent_workspace->variable_group(group_name);
+ GMM_ASSERT1(false, "Undefined variable group " << group_name);
+ }
+
+ const std::string&
+ ga_workspace::first_variable_of_group(const std::string &name) const {
+ const std::vector<std::string> &t = variable_group(name);
+ GMM_ASSERT1(t.size(), "Variable group " << name << " is empty");
+ return t[0];
+ }
+
+ bool ga_workspace::is_constant(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) return !(it->second.is_variable);
+ if (variable_group_exists(name))
+ return is_constant(first_variable_of_group(name));
+ if (md && md->variable_exists(name)) {
+ if (enable_all_md_variables) return md->is_true_data(name);
+ return md->is_data(name);
+ }
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->is_constant(name);
+ GMM_ASSERT1(false, "Undefined variable " << name);
+ }
+
+ bool ga_workspace::is_disabled_variable(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) return false;
+ if (variable_group_exists(name))
+ return is_disabled_variable(first_variable_of_group(name));
+ if (md && md->variable_exists(name)) {
+ if (enable_all_md_variables) return false;
+ return md->is_disabled_variable(name);
+ }
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->is_disabled_variable(name);
+ GMM_ASSERT1(false, "Undefined variable " << name);
+ }
+
+ const scalar_type &
+ ga_workspace::factor_of_variable(const std::string &name) const {
+ static const scalar_type one(1);
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) return one;
+ if (variable_group_exists(name))
+ return one;
+ if (md && md->variable_exists(name)) return md->factor_of_variable(name);
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->factor_of_variable(name);
+ GMM_ASSERT1(false, "Undefined variable " << name);
+ }
+
+ const gmm::sub_interval &
+ ga_workspace::interval_of_disabled_variable(const std::string &name) const {
+ std::map<std::string, gmm::sub_interval>::const_iterator
+ it1 = int_disabled_variables.find(name);
+ if (it1 != int_disabled_variables.end()) return it1->second;
+ if (md->is_affine_dependent_variable(name))
+ return interval_of_disabled_variable(md->org_variable(name));
+
+ size_type first = md->nb_dof();
+ for (const std::pair<std::string, gmm::sub_interval> &p
+ : int_disabled_variables)
+ first = std::max(first, p.second.last());
+
+ int_disabled_variables[name]
+ = gmm::sub_interval(first, gmm::vect_size(value(name)));
+ return int_disabled_variables[name];
+ }
+
+ const gmm::sub_interval &
+ ga_workspace::interval_of_variable(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) return it->second.I;
+ if (md && md->variable_exists(name)) {
+ if (enable_all_md_variables && md->is_disabled_variable(name))
+ return interval_of_disabled_variable(name);
+ return md->interval_of_variable(name);
+ }
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->interval_of_variable(name);
+ GMM_ASSERT1(false, "Undefined variable " << name);
+ }
+
+ const mesh_fem *
+ ga_workspace::associated_mf(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end())
+ return it->second.is_fem_dofs ? it->second.mf : 0;
+ if (md && md->variable_exists(name))
+ return md->pmesh_fem_of_variable(name);
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->associated_mf(name);
+ if (variable_group_exists(name))
+ return associated_mf(first_variable_of_group(name));
+ GMM_ASSERT1(false, "Undefined variable or group " << name);
+ }
+
+ const im_data *
+ ga_workspace::associated_im_data(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) return it->second.imd;
+ if (md && md->variable_exists(name))
+ return md->pim_data_of_variable(name);
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->associated_im_data(name);
+ if (variable_group_exists(name)) return 0;
+ GMM_ASSERT1(false, "Undefined variable " << name);
+ }
+
+ size_type ga_workspace::qdim(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) {
+ const mesh_fem *mf = it->second.is_fem_dofs ? it->second.mf : 0;
+ const im_data *imd = it->second.imd;
+ size_type n = it->second.qdim();
+ if (mf) {
+ return n * mf->get_qdim();
+ } else if (imd) {
+ return n * imd->tensor_size().total_size();
+ }
+ return n;
+ }
+ if (md && md->variable_exists(name))
+ return md->qdim_of_variable(name);
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->qdim(name);
+ if (variable_group_exists(name))
+ return qdim(first_variable_of_group(name));
+ GMM_ASSERT1(false, "Undefined variable or group " << name);
+ }
+
+ bgeot::multi_index
+ ga_workspace::qdims(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end()) {
+ const mesh_fem *mf = it->second.is_fem_dofs ? it->second.mf : 0;
+ const im_data *imd = it->second.imd;
+ size_type n = it->second.qdim();
+ if (mf) {
+ bgeot::multi_index mi = mf->get_qdims();
+ if (n > 1 || it->second.qdims.size() > 1) {
+ size_type i = 0;
+ if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
+ for (; i < it->second.qdims.size(); ++i)
+ mi.push_back(it->second.qdims[i]);
+ }
+ return mi;
+ } else if (imd) {
+ bgeot::multi_index mi = imd->tensor_size();
+ size_type q = n / imd->nb_filtered_index();
+ GMM_ASSERT1(q % imd->nb_tensor_elem() == 0,
+ "Invalid mesh im data vector");
+ if (n > 1 || it->second.qdims.size() > 1) {
+ size_type i = 0;
+ if (mi.back() == 1) { mi.back() *= it->second.qdims[0]; ++i; }
+ for (; i < it->second.qdims.size(); ++i)
+ mi.push_back(it->second.qdims[i]);
+ }
+ return mi;
+ }
+ return it->second.qdims;
+ }
+ if (md && md->variable_exists(name))
+ return md->qdims_of_variable(name);
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->qdims(name);
+ if (variable_group_exists(name))
+ return qdims(first_variable_of_group(name));
+ GMM_ASSERT1(false, "Undefined variable or group " << name);
+ }
+
+ const model_real_plain_vector &
+ ga_workspace::value(const std::string &name) const {
+ VAR_SET::const_iterator it = variables.find(name);
+ if (it != variables.end())
+ return *(it->second.V);
+ if (md && md->variable_exists(name))
+ return md->real_variable(name);
+ if (parent_workspace && parent_workspace->variable_exists(name))
+ return parent_workspace->value(name);
+ if (variable_group_exists(name))
+ return value(first_variable_of_group(name));
+ GMM_ASSERT1(false, "Undefined variable or group " << name);
+ }
+
+ scalar_type ga_workspace::get_time_step() const {
+ if (md) return md->get_time_step();
+ if (parent_workspace) return parent_workspace->get_time_step();
+ GMM_ASSERT1(false, "No time step defined here");
+ }
+
+
+ // Macros
+ bool ga_workspace::macro_exists(const std::string &name) const {
+ if (macros.find(name) != macros.end()) return true;
+ if (md && md->macro_exists(name)) return true;
+ if (parent_workspace &&
+ parent_workspace->macro_exists(name)) return true;
+ return false;
+ }
+
+ const std::string&
+ ga_workspace::get_macro(const std::string &name) const {
+ std::map<std::string, std::string>::const_iterator it=macros.find(name);
+ if (it != macros.end()) return it->second;
+ if (md && md->macro_exists(name)) return md->get_macro(name);
+ if (parent_workspace &&
+ parent_workspace->macro_exists(name))
+ return parent_workspace->get_macro(name);
+ GMM_ASSERT1(false, "Undefined macro");
+ }
+
+ ga_tree &
+ ga_workspace::macro_tree(const std::string &name, size_type meshdim,
+ size_type ref_elt_dim, bool ignore_X) const {
+ GMM_ASSERT1(macro_exists(name), "Undefined macro");
+ auto it = macro_trees.find(name);
+ bool to_be_analyzed = false;
+ m_tree *mt = 0;
+
+ if (it == macro_trees.end()) {
+ mt = &(macro_trees[name]);
+ to_be_analyzed = true;
+ } else {
+ mt = &(it->second);
+ GMM_ASSERT1(mt->ptree, "Recursive definition of macro " << name);
+ if (mt->meshdim != meshdim || mt->ignore_X != ignore_X) {
+ to_be_analyzed = true;
+ delete mt->ptree; mt->ptree = 0;
+ }
+ }
+ if (to_be_analyzed) {
+ ga_tree tree;
+ ga_read_string(get_macro(name), tree);
+ ga_semantic_analysis(get_macro(name), tree, *this, meshdim, ref_elt_dim,
+ false, ignore_X, 3);
+ GMM_ASSERT1(tree.root, "Invalid macro");
+ mt->ptree = new ga_tree(tree);
+ mt->meshdim = meshdim;
+ mt->ignore_X = ignore_X;
+ }
+ return *(mt->ptree);
+ }
+
+ void ga_workspace::add_interpolate_transformation
+ (const std::string &name, pinterpolate_transformation ptrans) {
+ if (transformations.find(name) != transformations.end())
+ GMM_ASSERT1(name.compare("neighbour_elt"), "neighbour_elt is a "
+ "reserved interpolate transformation name");
+ transformations[name] = ptrans;
+ }
+
+ bool ga_workspace::interpolate_transformation_exists
+ (const std::string &name) const {
+ return (md && md->interpolate_transformation_exists(name)) ||
+ (parent_workspace &&
+ parent_workspace->interpolate_transformation_exists(name)) ||
+ (transformations.find(name) != transformations.end());
+ }
+
+ pinterpolate_transformation
+ ga_workspace::interpolate_transformation(const std::string &name) const {
+ std::map<std::string, pinterpolate_transformation>::const_iterator
+ it = transformations.find(name);
+ if (it != transformations.end()) return it->second;
+ if (md && md->interpolate_transformation_exists(name))
+ return md->interpolate_transformation(name);
+ if (parent_workspace &&
+ parent_workspace->interpolate_transformation_exists(name))
+ return parent_workspace->interpolate_transformation(name);
+ GMM_ASSERT1(false, "Inexistent transformation " << name);
+ }
+
+ bool ga_workspace::elementary_transformation_exists
+ (const std::string &name) const {
+ return (md && md->elementary_transformation_exists(name)) ||
+ (parent_workspace &&
+ parent_workspace->elementary_transformation_exists(name)) ||
+ (elem_transformations.find(name) != elem_transformations.end());
+ }
+
+ pelementary_transformation
+ ga_workspace::elementary_transformation(const std::string &name) const {
+ std::map<std::string, pelementary_transformation>::const_iterator
+ it = elem_transformations.find(name);
+ if (it != elem_transformations.end()) return it->second;
+ if (md && md->elementary_transformation_exists(name))
+ return md->elementary_transformation(name);
+ if (parent_workspace &&
+ parent_workspace->elementary_transformation_exists(name))
+ return parent_workspace->elementary_transformation(name);
+ GMM_ASSERT1(false, "Inexistent elementary transformation " << name);
+ }
+
+ const mesh_region &
+ ga_workspace::register_region(const mesh &m, const mesh_region ®ion) {
+ if (&m == &dummy_mesh())
+ return dummy_mesh_region();
+
+ std::list<mesh_region> &lmr = registred_mesh_regions[&m];
+ for (const mesh_region &rg : lmr)
+ if (rg.compare(m, region, m)) return rg;
+ lmr.push_back(region);
+ return lmr.back();
+ }
+
+
+ void ga_workspace::add_tree(ga_tree &tree, const mesh &m,
+ const mesh_im &mim, const mesh_region &rg,
+ const std::string &expr,
+ size_type add_derivative_order,
+ bool function_expr, size_type for_interpolation,
+ const std::string varname_interpolation) {
+ if (tree.root) {
+
+ // Eliminate the term if it corresponds to disabled variables
+ if ((tree.root->test_function_type >= 1 &&
+ is_disabled_variable(tree.root->name_test1)) ||
+ (tree.root->test_function_type >= 2 &&
+ is_disabled_variable(tree.root->name_test2))) {
+ // cout<<"disabling term "; ga_print_node(tree.root, cout);
cout<<endl;
+ return;
+ }
+ // cout << "add tree with tests functions of " << tree.root->name_test1
+ // << " and " << tree.root->name_test2 << endl;
+ // ga_print_node(tree.root, cout); cout << endl;
+ bool remain = true;
+ size_type order = 0, ind_tree = 0;
+
+ if (for_interpolation)
+ order = size_type(-1) - add_derivative_order;
+ else {
+ switch(tree.root->test_function_type) {
+ case 0: order = 0; break;
+ case 1: order = 1; break;
+ case 3: order = 2; break;
+ default: GMM_ASSERT1(false, "Inconsistent term "
+ << tree.root->test_function_type);
+ }
+ }
+
+ bool found = false;
+ for (size_type i = 0; i < trees.size(); ++i) {
+ if (trees[i].mim == &mim && trees[i].m == &m &&
+ trees[i].order == order &&
+ trees[i].name_test1.compare(tree.root->name_test1) == 0 &&
+ trees[i].interpolate_name_test1.compare
+ (tree.root->interpolate_name_test1) == 0 &&
+ trees[i].name_test2.compare(tree.root->name_test2) == 0 &&
+ trees[i].interpolate_name_test2.compare
+ (tree.root->interpolate_name_test2) == 0 &&
+ trees[i].rg == &rg && trees[i].interpolation == for_interpolation
&&
+ trees[i].varname_interpolation.compare(varname_interpolation)==0) {
+ ga_tree &ftree = *(trees[i].ptree);
+
+ ftree.insert_node(ftree.root, GA_NODE_OP);
+ ftree.root->op_type = GA_PLUS;
+ ftree.root->children.resize(2, nullptr);
+ ftree.copy_node(tree.root, ftree.root, ftree.root->children[1]);
+ ga_semantic_analysis("", ftree, *this, m.dim(),
+ ref_elt_dim_of_mesh(m), false, function_expr);
+ found = true;
+ break;
+ }
+ }
+
+ if (!found) {
+ ind_tree = trees.size(); remain = false;
+ trees.push_back(tree_description());
+ trees.back().mim = &mim; trees.back().m = &m;
+ trees.back().rg = &rg;
+ trees.back().ptree = new ga_tree;
+ trees.back().ptree->swap(tree);
+ pga_tree_node root = trees.back().ptree->root;
+ trees.back().name_test1 = root->name_test1;
+ trees.back().name_test2 = root->name_test2;
+ trees.back().interpolate_name_test1 = root->interpolate_name_test1;
+ trees.back().interpolate_name_test2 = root->interpolate_name_test2;
+ trees.back().order = order;
+ trees.back().interpolation = for_interpolation;
+ trees.back().varname_interpolation = varname_interpolation;
+ }
+
+ if (for_interpolation == 0 && order < add_derivative_order) {
+ std::set<var_trans_pair> expr_variables;
+ ga_extract_variables((remain ? tree : *(trees[ind_tree].ptree)).root,
+ *this, m, expr_variables, true);
+ for (const var_trans_pair &var : expr_variables) {
+ if (!(is_constant(var.varname))) {
+ ga_tree dtree = (remain ? tree : *(trees[ind_tree].ptree));
+ // cout << "Derivation with respect to " << var.first << " : "
+ // << var.second << " of " << ga_tree_to_string(dtree) << endl;
+ GA_TIC;
+ ga_derivative(dtree, *this, m, var.varname, var.transname,
1+order);
+ // cout << "Result : " << ga_tree_to_string(dtree) << endl;
+ GA_TOCTIC("Derivative time");
+ ga_semantic_analysis(expr, dtree, *this, m.dim(),
+ ref_elt_dim_of_mesh(m), false, function_expr);
+ GA_TOCTIC("Analysis after Derivative time");
+ // cout << "after analysis " << ga_tree_to_string(dtree) << endl;
+ add_tree(dtree, m, mim, rg, expr, add_derivative_order,
+ function_expr, for_interpolation, varname_interpolation);
+ }
+ }
+ }
+ }
+ }
+
+ ga_workspace::m_tree::~m_tree() { if (ptree) delete ptree; }
+ ga_workspace::m_tree::m_tree(const m_tree& o)
+ : ptree(o.ptree), meshdim(o.meshdim), ignore_X(o.ignore_X)
+ { if (o.ptree) ptree = new ga_tree(*(o.ptree)); }
+ ga_workspace::m_tree &ga_workspace::m_tree::operator =(const m_tree& o) {
+ if (ptree) delete ptree;
+ ptree = o.ptree; meshdim = o.meshdim; ignore_X = o.ignore_X;
+ if (o.ptree) ptree = new ga_tree(*(o.ptree));
+ return *this;
+ }
+
+ size_type ga_workspace::add_expression(const std::string &expr,
+ const mesh_im &mim,
+ const mesh_region &rg_,
+ size_type add_derivative_order) {
+ const mesh_region &rg = register_region(mim.linked_mesh(), rg_);
+ // cout << "adding expression " << expr << endl;
+ GA_TIC;
+ size_type max_order = 0;
+ std::vector<ga_tree> ltrees(1);
+ ga_read_string(expr, ltrees[0]);
+ // cout << "read : " << ga_tree_to_string(ltrees[0]) << endl;
+ ga_semantic_analysis(expr, ltrees[0], *this, mim.linked_mesh().dim(),
+ ref_elt_dim_of_mesh(mim.linked_mesh()),
+ false, false, 1);
+ // cout << "analysed : " << ga_tree_to_string(ltrees[0]) << endl;
+ GA_TOC("First analysis time");
+ if (ltrees[0].root) {
+ if (test1.size() > 1 || test2.size() > 1) {
+ size_type ntest2 = std::max(size_type(1), test2.size());
+ size_type nb_ltrees = test1.size()*ntest2;
+ ltrees.resize(nb_ltrees);
+ for (size_type i = 1; i < nb_ltrees; ++i) ltrees[i] = ltrees[0];
+ std::set<var_trans_pair>::iterator it1 = test1.begin();
+ for (size_type i = 0; i < test1.size(); ++i, ++it1) {
+ std::set<var_trans_pair>::iterator it2 = test2.begin();
+ for (size_type j = 0; j < ntest2; ++j) {
+ selected_test1 = *it1;
+ if (test2.size()) selected_test2 = *it2++;
+ // cout << "analysis with " << selected_test1.first << endl;
+ ga_semantic_analysis(expr, ltrees[i*ntest2+j], *this,
+ mim.linked_mesh().dim(),
+ ref_elt_dim_of_mesh(mim.linked_mesh()),
+ false, false, 2);
+ // cout <<"split: "<< ga_tree_to_string(ltrees[i*ntest2+j]) <<
endl;
+ }
+ }
+ }
+
+ for (size_type i = 0; i < ltrees.size(); ++i) {
+ if (ltrees[i].root) {
+ // cout << "adding tree " << ga_tree_to_string(ltrees[i]) << endl;
+ max_order = std::max(ltrees[i].root->nb_test_functions(), max_order);
+ add_tree(ltrees[i], mim.linked_mesh(), mim, rg, expr,
+ add_derivative_order, true, 0, "");
+ }
+ }
+ }
+ GA_TOC("Time for add expression");
+ return max_order;
+ }
+
+ void ga_workspace::add_function_expression(const std::string &expr) {
+ ga_tree tree;
+ ga_read_string(expr, tree);
+ ga_semantic_analysis(expr, tree, *this, 1, 1, false, true);
+ if (tree.root) {
+ // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
+ // "Invalid function expression");
+ add_tree(tree, dummy_mesh(), dummy_mesh_im(), dummy_mesh_region(),
+ expr, 0, true, 0, "");
+ }
+ }
+
+ void ga_workspace::add_interpolation_expression(const std::string &expr,
+ const mesh &m,
+ const mesh_region &rg_) {
+ const mesh_region &rg = register_region(m, rg_);
+ ga_tree tree;
+ ga_read_string(expr, tree);
+ ga_semantic_analysis(expr, tree, *this, m.dim(), ref_elt_dim_of_mesh(m),
+ false, false);
+ if (tree.root) {
+ // GMM_ASSERT1(tree.root->nb_test_functions() == 0,
+ // "Invalid expression containing test functions");
+ add_tree(tree, m, dummy_mesh_im(), rg, expr, 0, false, 1, "");
+ }
+ }
+
+ void ga_workspace::add_interpolation_expression(const std::string &expr,
+ const mesh_im &mim,
+ const mesh_region &rg_) {
+ const mesh &m = mim.linked_mesh();
+ const mesh_region &rg = register_region(m, rg_);
+ ga_tree tree;
+ ga_read_string(expr, tree);
+ ga_semantic_analysis(expr, tree, *this, m.dim(), ref_elt_dim_of_mesh(m),
+ false, false);
+ if (tree.root) {
+ GMM_ASSERT1(tree.root->nb_test_functions() == 0,
+ "Invalid expression containing test functions");
+ add_tree(tree, m, mim, rg, expr, 0, false, 1, "");
+ }
+ }
+
+ void ga_workspace::add_assignment_expression
+ (const std::string &varname, const std::string &expr, const mesh_region &rg_,
+ size_type order, bool before) {
+ const im_data *imd = associated_im_data(varname);
+ GMM_ASSERT1(imd != 0, "Only applicable to im_data");
+ const mesh_im &mim = imd->linked_mesh_im();
+ const mesh &m = mim.linked_mesh();
+ const mesh_region &rg = register_region(m, rg_);
+ ga_tree tree;
+ ga_read_string(expr, tree);
+ ga_semantic_analysis(expr, tree, *this, m.dim(), ref_elt_dim_of_mesh(m),
+ false, false);
+ if (tree.root) {
+ GMM_ASSERT1(tree.root->nb_test_functions() == 0,
+ "Invalid expression containing test functions");
+ add_tree(tree, m, mim, rg, expr, order+1, false, (before ? 1 : 2),
+ varname);
+ }
+ }
+
+ size_type ga_workspace::nb_trees() const { return trees.size(); }
+
+ ga_workspace::tree_description &ga_workspace::tree_info(size_type i)
+ { return trees[i]; }
+
+ bool ga_workspace::used_variables(std::vector<std::string> &vl,
+ std::vector<std::string> &vl_test1,
+ std::vector<std::string> &vl_test2,
+ std::vector<std::string> &dl,
+ size_type order) {
+ bool islin = true;
+ std::set<var_trans_pair> vll, dll;
+ for (size_type i = 0; i < vl.size(); ++i)
+ vll.insert(var_trans_pair(vl[i], ""));
+ for (size_type i = 0; i < dl.size(); ++i)
+ dll.insert(var_trans_pair(dl[i], ""));
+
+ for (size_type i = 0; i < trees.size(); ++i) {
+ ga_workspace::tree_description &td = trees[i];
+ std::set<var_trans_pair> dllaux;
+ bool fv = ga_extract_variables(td.ptree->root, *this, *(td.m),
+ dllaux, false);
+
+ if (td.order == order) {
+ for (std::set<var_trans_pair>::iterator it = dllaux.begin();
+ it!=dllaux.end(); ++it)
+ dll.insert(*it);
+ }
+ switch (td.order) {
+ case 0: break;
+ case 1:
+ if (td.order == order) {
+ if (variable_group_exists(td.name_test1)) {
+ for (const std::string &t : variable_group(td.name_test1))
+ vll.insert(var_trans_pair(t, td.interpolate_name_test1));
+ } else {
+ vll.insert(var_trans_pair(td.name_test1,
+ td.interpolate_name_test1));
+ }
+ bool found = false;
+ for (const std::string &t : vl_test1)
+ if (td.name_test1.compare(t) == 0)
+ found = true;
+ if (!found)
+ vl_test1.push_back(td.name_test1);
+ }
+ break;
+ case 2:
+ if (td.order == order) {
+ if (variable_group_exists(td.name_test1)) {
+ for (const std::string &t : variable_group(td.name_test1))
+ vll.insert(var_trans_pair(t, td.interpolate_name_test1));
+ } else {
+ vll.insert(var_trans_pair(td.name_test1,
+ td.interpolate_name_test1));
+ }
+ if (variable_group_exists(td.name_test2)) {
+ for (const std::string &t : variable_group(td.name_test2))
+ vll.insert(var_trans_pair(t, td.interpolate_name_test2));
+ } else {
+ vll.insert(var_trans_pair(td.name_test2,
+ td.interpolate_name_test2));
+ }
+ bool found = false;
+ for (size_type j = 0; j < vl_test1.size(); ++j)
+ if ((td.name_test1.compare(vl_test1[j]) == 0) &&
+ (td.name_test2.compare(vl_test2[j]) == 0))
+ found = true;
+ if (!found) {
+ vl_test1.push_back(td.name_test1);
+ vl_test2.push_back(td.name_test2);
+ }
+ }
+ if (fv) islin = false;
+ break;
+ }
+ }
+ vl.clear();
+ for (const auto &var : vll)
+ if (vl.size() == 0 || var.varname.compare(vl.back()))
+ vl.push_back(var.varname);
+ dl.clear();
+ for (const auto &var : dll)
+ if (dl.size() == 0 || var.varname.compare(dl.back()))
+ dl.push_back(var.varname);
+
+ return islin;
+ }
+
+ void ga_workspace::define_variable_group(const std::string &group_name,
+ const std::vector<std::string> &nl)
{
+ GMM_ASSERT1(!(variable_exists(group_name)), "The name of a group of "
+ "variables cannot be the same as a variable name");
+
+ std::set<const mesh *> ms;
+ bool is_data_ = false;
+ for (size_type i = 0; i < nl.size(); ++i) {
+ if (i == 0)
+ is_data_ = is_constant(nl[i]);
+ else {
+ GMM_ASSERT1(is_data_ == is_constant(nl[i]),
+ "It is not possible to mix variables and data in a group");
+ }
+ GMM_ASSERT1(variable_exists(nl[i]),
+ "All variables in a group have to exist in the model");
+ const mesh_fem *mf = associated_mf(nl[i]);
+ GMM_ASSERT1(mf, "Variables in a group should be fem variables");
+ GMM_ASSERT1(ms.find(&(mf->linked_mesh())) == ms.end(),
+ "Two variables in a group cannot share the same mesh");
+ ms.insert(&(mf->linked_mesh()));
+ }
+ variable_groups[group_name] = nl;
+ }
+
+
+ const std::string &ga_workspace::variable_in_group
+ (const std::string &group_name, const mesh &m) const {
+ if (variable_group_exists(group_name)) {
+ for (const std::string &t : variable_group(group_name))
+ if (&(associated_mf(t)->linked_mesh()) == &m)
+ return t;
+ GMM_ASSERT1(false, "No variable in this group for the given mesh");
+ } else
+ return group_name;
+ }
+
+
+ void ga_workspace::assembly(size_type order) {
+ size_type ndof;
+ const ga_workspace *w = this;
+ while (w->parent_workspace) w = w->parent_workspace;
+ if (w->md) ndof = w->md->nb_dof(); // To eventually call actualize_sizes()
+
+ GA_TIC;
+ ga_instruction_set gis;
+ ga_compile(*this, gis, order);
+ ndof = gis.nb_dof;
+ size_type max_dof = gis.max_dof;
+ GA_TOCTIC("Compile time");
+
+ if (order == 2) {
+ if (K.use_count()) {
+ gmm::clear(*K);
+ gmm::resize(*K, max_dof, max_dof);
+ }
+ gmm::clear(unreduced_K);
+ gmm::resize(unreduced_K, ndof, ndof);
+ }
+ if (order == 1) {
+ if (V.use_count()) {
+ gmm::clear(*V);
+ gmm::resize(*V, max_dof);
+ }
+ gmm::clear(unreduced_V);
+ gmm::resize(unreduced_V, ndof);
+ }
+ E = 0;
+ GA_TOCTIC("Init time");
+
+ ga_exec(gis, *this);
+ GA_TOCTIC("Exec time");
+
+ if (order == 1) {
+ MPI_SUM_VECTOR(assembled_vector());
+ MPI_SUM_VECTOR(unreduced_V);
+ } else if (order == 0) {
+ assembled_potential() = MPI_SUM_SCALAR(assembled_potential());
+ }
+
+ // Deal with reduced fems.
+ if (order > 0) {
+ std::set<std::string> vars_vec_done;
+ std::set<std::pair<std::string, std::string> > vars_mat_done;
+ for (ga_tree &tree : gis.trees) {
+ if (tree.root) {
+ if (order == 1) {
+ const std::string &name = tree.root->name_test1;
+ const std::vector<std::string> vnames_(1,name);
+ const std::vector<std::string> &vnames
+ = variable_group_exists(name) ? variable_group(name)
+ : vnames_;
+ for (const std::string &vname : vnames) {
+ const mesh_fem *mf = associated_mf(vname);
+ if (mf && mf->is_reduced() &&
+ vars_vec_done.find(vname) == vars_vec_done.end()) {
+ gmm::mult_add(gmm::transposed(mf->extension_matrix()),
+ gmm::sub_vector(unreduced_V,
+ gis.var_intervals[vname]),
+ gmm::sub_vector(*V,
+ interval_of_variable(vname)));
+ vars_vec_done.insert(vname);
+ }
+ }
+ } else {
+ std::string &name1 = tree.root->name_test1;
+ std::string &name2 = tree.root->name_test2;
+ const std::vector<std::string> vnames1_(1,name1),
+ vnames2_(2,name2);
+ const std::vector<std::string> &vnames1
+ = variable_group_exists(name1) ? variable_group(name1)
+ : vnames1_;
+ const std::vector<std::string> &vnames2
+ = variable_group_exists(name2) ? variable_group(name2)
+ : vnames2_;
+ for (const std::string &vname1 : vnames1) {
+ for (const std::string &vname2 : vnames2) {
+ const mesh_fem *mf1 = associated_mf(vname1);
+ const mesh_fem *mf2 = associated_mf(vname2);
+ if (((mf1 && mf1->is_reduced())
+ || (mf2 && mf2->is_reduced()))) {
+ std::pair<std::string, std::string> p(vname1, vname2);
+ if (vars_mat_done.find(p) == vars_mat_done.end()) {
+ gmm::sub_interval uI1 = gis.var_intervals[vname1];
+ gmm::sub_interval uI2 = gis.var_intervals[vname2];
+ gmm::sub_interval I1 = interval_of_variable(vname1);
+ gmm::sub_interval I2 = interval_of_variable(vname2);
+ if ((mf1 && mf1->is_reduced()) &&
+ (mf2 && mf2->is_reduced())) {
+ model_real_sparse_matrix aux(I1.size(), uI2.size());
+ model_real_row_sparse_matrix M(I1.size(), I2.size());
+ gmm::mult(gmm::transposed(mf1->extension_matrix()),
+ gmm::sub_matrix(unreduced_K, uI1, uI2), aux);
+ gmm::mult(aux, mf2->extension_matrix(), M);
+ gmm::add(M, gmm::sub_matrix(*K, I1, I2));
+ } else if (mf1 && mf1->is_reduced()) {
+ model_real_sparse_matrix M(I1.size(), I2.size());
+ gmm::mult(gmm::transposed(mf1->extension_matrix()),
+ gmm::sub_matrix(unreduced_K, uI1, uI2), M);
+ gmm::add(M, gmm::sub_matrix(*K, I1, I2));
+ } else {
+ model_real_row_sparse_matrix M(I1.size(), I2.size());
+ gmm::mult(gmm::sub_matrix(unreduced_K, uI1, uI2),
+ mf2->extension_matrix(), M);
+ gmm::add(M, gmm::sub_matrix(*K, I1, I2));
+ }
+ vars_mat_done.insert(p);
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ void ga_workspace::clear_expressions() {
+ trees.clear();
+ macro_trees.clear();
+ }
+
+ void ga_workspace::print(std::ostream &str) {
+ for (size_type i = 0; i < trees.size(); ++i)
+ if (trees[i].ptree->root) {
+ cout << "Expression tree " << i << " of order " <<
+ trees[i].ptree->root->nb_test_functions() << " :" << endl;
+ ga_print_node(trees[i].ptree->root, str);
+ cout << endl;
+ }
+ }
+
+ void ga_workspace::tree_description::copy(const tree_description& td) {
+ order = td.order;
+ interpolation = td.interpolation;
+ varname_interpolation = td.varname_interpolation;
+ name_test1 = td.name_test1;
+ name_test2 = td.name_test2;
+ interpolate_name_test1 = td.interpolate_name_test1;
+ interpolate_name_test2 = td.interpolate_name_test2;
+ mim = td.mim;
+ m = td.m;
+ rg = td.rg;
+ ptree = 0;
+ if (td.ptree) ptree = new ga_tree(*(td.ptree));
+ }
+
+ ga_workspace::tree_description &ga_workspace::tree_description::operator =
+ (const ga_workspace::tree_description& td)
+ { if (ptree) delete ptree; ptree = 0; copy(td); return *this; }
+ ga_workspace::tree_description::~tree_description()
+ { if (ptree) delete ptree; ptree = 0; }
+
+
+
+ //=========================================================================
+ // Extract the constant term of degree 1 expressions
+ //=========================================================================
+
+ std::string ga_workspace::extract_constant_term(const mesh &m) {
+ std::string constant_term;
+ for (size_type i = 0; i < trees.size(); ++i) {
+ ga_workspace::tree_description &td = trees[i];
+
+ if (td.order == 1) {
+ ga_tree local_tree = *(td.ptree);
+ if (local_tree.root)
+ ga_node_extract_constant_term(local_tree, local_tree.root, *this, m);
+ if (local_tree.root)
+ ga_semantic_analysis("", local_tree, *this, m.dim(),
+ ref_elt_dim_of_mesh(m), false, false);
+ if (local_tree.root && local_tree.root->node_type != GA_NODE_ZERO) {
+ constant_term += "-("+ga_tree_to_string(local_tree)+")";
+ }
+ }
+ }
+ return constant_term;
+ }
+
+ //=========================================================================
+ // Extract the order zero term
+ //=========================================================================
+
+ std::string ga_workspace::extract_order0_term() {
+ std::string term;
+ for (size_type i = 0; i < trees.size(); ++i) {
+ ga_workspace::tree_description &td = trees[i];
+ if (td.order == 0) {
+ ga_tree &local_tree = *(td.ptree);
+ if (term.size())
+ term += "+("+ga_tree_to_string(local_tree)+")";
+ else
+ term = "("+ga_tree_to_string(local_tree)+")";
+ }
+ }
+ return term;
+ }
+
+
+ //=========================================================================
+ // Extract the order one term corresponding to a certain test function
+ //=========================================================================
+
+ std::string ga_workspace::extract_order1_term(const std::string &varname) {
+ std::string term;
+ for (size_type i = 0; i < trees.size(); ++i) {
+ ga_workspace::tree_description &td = trees[i];
+ if (td.order == 1 && td.name_test1.compare(varname) == 0) {
+ ga_tree &local_tree = *(td.ptree);
+ if (term.size())
+ term += "+("+ga_tree_to_string(local_tree)+")";
+ else
+ term = "("+ga_tree_to_string(local_tree)+")";
+ }
+ }
+ return term;
+ }
+
+ //=========================================================================
+ // Extract Neumann terms
+ //=========================================================================
+
+ std::string ga_workspace::extract_Neumann_term(const std::string &varname) {
+ std::string result;
+ for (size_type i = 0; i < trees.size(); ++i) {
+ ga_workspace::tree_description &td = trees[i];
+ if (td.order == 1 && td.name_test1.compare(varname) == 0) {
+ ga_tree &local_tree = *(td.ptree);
+ if (local_tree.root)
+ ga_extract_Neumann_term(local_tree, varname, *this,
+ local_tree.root, result);
+ }
+ }
+ return result;
+ }
+
+} /* end of namespace */