[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Sun, 3 Mar 2019 17:17:20 -0500 (EST) |
branch: master
commit 2f7396a46dadad266ba882623d300d1a19e274dc
Author: Konstantinos Poulios <address@hidden>
Date: Sun Mar 3 23:17:02 2019 +0100
Code cleanup and stricter constness
---
src/getfem/getfem_generic_assembly.h | 42 ++++-----
.../getfem_generic_assembly_compile_and_exec.h | 25 ++----
src/getfem_generic_assembly_compile_and_exec.cc | 100 ++++++++-------------
src/getfem_generic_assembly_interpolation.cc | 19 +++-
src/getfem_generic_assembly_workspace.cc | 52 ++++++-----
5 files changed, 111 insertions(+), 127 deletions(-)
diff --git a/src/getfem/getfem_generic_assembly.h
b/src/getfem/getfem_generic_assembly.h
index 5b91109..4312ec8 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -270,13 +270,14 @@ namespace getfem {
struct var_description {
- bool is_variable;
- bool is_fem_dofs;
+ const bool is_variable;
+ const bool is_fem_dofs;
const mesh_fem *mf;
+ const im_data *imd;
gmm::sub_interval I;
const model_real_plain_vector *V;
- const im_data *imd;
- bgeot::multi_index qdims; // For data having a qdim != of the fem
+ bgeot::multi_index qdims; // For data having a qdim different than
+ // the qdim of the fem or im_data
// (dim per dof for dof data)
// and for constant variables.
@@ -286,28 +287,29 @@ namespace getfem {
return q;
}
- var_description(bool is_var, bool is_fem,
- const mesh_fem *mmf, gmm::sub_interval I_,
- const model_real_plain_vector *v, const im_data *imd_,
+ var_description(bool is_var, const mesh_fem *mf_, const im_data *imd_,
+ gmm::sub_interval I_, const model_real_plain_vector *v,
size_type Q)
- : is_variable(is_var), is_fem_dofs(is_fem), mf(mmf), I(I_), V(v),
- imd(imd_), qdims(1) {
+ : is_variable(is_var), is_fem_dofs(mf_ != 0), mf(mf_), imd(imd_),
+ I(I_), V(v), qdims(1)
+ {
GMM_ASSERT1(Q > 0, "Bad dimension");
qdims[0] = Q;
}
- var_description() : is_variable(false), is_fem_dofs(false),
- mf(0), V(0), imd(0), qdims(1) { qdims[0] = 1; }
};
public:
+ enum operation_type {ASSEMBLY,
+ PRE_ASSIGNMENT,
+ POST_ASSIGNMENT};
+
struct tree_description { // CAUTION: Specific copy constructor
- size_type order; // 0: potential, 1: weak form, 2: tangent operator
- // -1 : interpolation/ assignment all order,
- // -2 : assignment on potential, -3 : assignment on weak form
- // -3 : assignment on tangent operator
- size_type interpolation; // O : assembly, 1 : interpolate before assembly
- // 2 : interpolate after assembly.
+ size_type order; // 0 : potential
+ // 1 : residual
+ // 2 : tangent operator
+ // -1 : any
+ operation_type operation;
std::string varname_interpolation; // Where to interpolate
std::string name_test1, name_test2;
std::string interpolate_name_test1, interpolate_name_test2;
@@ -317,7 +319,7 @@ namespace getfem {
const mesh_region *rg;
ga_tree *ptree;
tree_description()
- : interpolation(0), varname_interpolation(""),
+ : operation(ASSEMBLY), varname_interpolation(""),
name_test1(""), name_test2(""),
interpolate_name_test1(""), interpolate_name_test2(""),
mim(0), m(0), rg(0), ptree(0) {}
@@ -366,8 +368,8 @@ namespace getfem {
void 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 scalar_expr, size_type for_interpolation,
- const std::string varname_interpolation);
+ bool scalar_expr, operation_type op_type=ASSEMBLY,
+ const std::string varname_interpolation="");
std::shared_ptr<model_real_sparse_matrix> K;
diff --git a/src/getfem/getfem_generic_assembly_compile_and_exec.h
b/src/getfem/getfem_generic_assembly_compile_and_exec.h
index 36588f5..888b2ce 100644
--- a/src/getfem/getfem_generic_assembly_compile_and_exec.h
+++ b/src/getfem/getfem_generic_assembly_compile_and_exec.h
@@ -69,9 +69,7 @@ namespace getfem {
};
typedef std::shared_ptr<ga_instruction> pga_instruction;
- typedef std::vector<pga_instruction> ga_instruction_list;
-
struct gauss_pt_corresp { // For neighbour interpolation transformation
bgeot::pgeometric_trans pgt1, pgt2;
papprox_integration pai;
@@ -185,25 +183,22 @@ namespace getfem {
std::map<std::string, elementary_trans_info> elementary_trans_infos;
secondary_domain_info secondary_domain_infos;
- // Instructions being executed at the first Gauss point after
- // a change of integration method only.
- ga_instruction_list begin_instructions;
- // Instructions executed once per element
- ga_instruction_list elt_instructions;
- // Instructions executed on each integration/interpolation point
- ga_instruction_list instructions;
+ std::vector<pga_instruction>
+ begin_instructions, // Instructions being executed at the first Gauss
+ // point after a change of integration method
only.
+ elt_instructions, // Instructions executed once per element
+ instructions; // Instructions executed on each
+ // integration/interpolation point
std::map<scalar_type, std::list<pga_tree_node> > node_list;
- region_mim_instructions(): m(0), im(0) {}
+ region_mim_instructions(): m(0), im(0) {}
};
std::list<ga_tree> trees; // The trees are stored mainly because they
// contain the intermediary tensors.
std::list<ga_tree> interpolation_trees;
- typedef std::map<region_mim, region_mim_instructions> instructions_set;
-
- instructions_set whole_instructions;
+ std::map<region_mim, region_mim_instructions> all_instructions;
ga_instruction_set() { max_dof = nb_dof = 0; need_elt_size = false; ipt=0;
}
};
@@ -220,10 +215,6 @@ namespace getfem {
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_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 6f96e86..cb7bd65 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -6592,17 +6592,17 @@ namespace getfem {
GMM_ASSERT1(!scalar || (root->tensor().size() == 1),
"The result of the given expression is not a scalar");
ga_instruction_set::region_mim rm(td.mim, td.rg, 0);
- gis.whole_instructions[rm].m = td.m;
+ gis.all_instructions[rm].m = td.m;
ga_if_hierarchy if_hierarchy;
- ga_compile_node(root, workspace, gis,
- gis.whole_instructions[rm],*(td.m),true,if_hierarchy);
+ ga_compile_node(root, workspace, gis, gis.all_instructions[rm],
+ *(td.m), true, if_hierarchy);
gis.coeff = scalar_type(1);
pga_instruction pgai;
workspace.assembled_tensor() = root->tensor();
pgai = std::make_shared<ga_instruction_add_to_coeff>
(workspace.assembled_tensor(), root->tensor(), gis.coeff);
- gis.whole_instructions[rm].instructions.push_back(std::move(pgai));
+ gis.all_instructions[rm].instructions.push_back(std::move(pgai));
}
}
}
@@ -6696,10 +6696,10 @@ namespace getfem {
void ga_compile_interpolation(ga_workspace &workspace,
ga_instruction_set &gis) {
gis.transformations.clear();
- gis.whole_instructions.clear();
+ gis.all_instructions.clear();
for (size_type i = 0; i < workspace.nb_trees(); ++i) {
const ga_workspace::tree_description &td = workspace.tree_info(i);
- if (td.interpolation > 0) {
+ if (td.operation != ga_workspace::ASSEMBLY) {
gis.trees.push_back(*(td.ptree));
// Semantic analysis mainly to evaluate fixed size variables and data
@@ -6711,8 +6711,7 @@ namespace getfem {
if (root) {
// Compile tree
ga_instruction_set::region_mim rm(td.mim, td.rg, 0);
- ga_instruction_set::region_mim_instructions &rmi
- = gis.whole_instructions[rm];
+ auto &rmi = gis.all_instructions[rm];
rmi.m = td.m;
rmi.im = td.mim;
// rmi.interpolate_infos.clear();
@@ -6733,40 +6732,37 @@ namespace getfem {
void ga_compile(ga_workspace &workspace,
ga_instruction_set &gis, size_type order) {
gis.transformations.clear();
- gis.whole_instructions.clear();
- for (size_type version : std::array<size_type, 3>{1, 0, 2}) {
+ gis.all_instructions.clear();
+ std::array<ga_workspace::operation_type,3>
+ phases{ga_workspace::PRE_ASSIGNMENT,
+ ga_workspace::ASSEMBLY,
+ ga_workspace::POST_ASSIGNMENT};
+ for (const auto &phase : phases) {
+
for (size_type i = 0; i < workspace.nb_trees(); ++i) {
ga_workspace::tree_description &td = workspace.tree_info(i);
-
- if ((version == td.interpolation) &&
- ((version == 0 && td.order == order) || // Assembly
- ((version > 0 && (td.order == size_type(-1) || // Assignment
- td.order == size_type(-2) - order))))) {
- ga_tree *added_tree = 0;
- if (td.interpolation) {
- gis.interpolation_trees.push_back(*(td.ptree));
- added_tree = &(gis.interpolation_trees.back());
- } else {
- gis.trees.push_back(*(td.ptree));
- added_tree = &(gis.trees.back());
- }
-
+ if (td.operation != phase)
+ continue; // skip this tree in this phase
+
+ if (td.order == order || td.order == size_type(-1)) {
+ std::list<ga_tree> &trees = (phase == ga_workspace::ASSEMBLY)
+ ? gis.trees
+ : gis.interpolation_trees;
+ trees.push_back(*(td.ptree));
// Semantic analysis mainly to evaluate fixed size variables and data
- ga_semantic_analysis(*added_tree, workspace,
- td.mim->linked_mesh(),
+ ga_semantic_analysis(trees.back(), workspace, td.mim->linked_mesh(),
ref_elt_dim_of_mesh(td.mim->linked_mesh()),
true, false);
- pga_tree_node root = added_tree->root;
+ pga_tree_node root = trees.back().root;
if (root) {
// Compile tree
// cout << "Will compile "; ga_print_node(root, cout); cout <<
endl;
psecondary_domain psd(0);
- if (added_tree->secondary_domain.size())
- psd = workspace.secondary_domain(added_tree->secondary_domain);
+ if (trees.back().secondary_domain.size())
+ psd = workspace.secondary_domain(trees.back().secondary_domain);
ga_instruction_set::region_mim rm(td.mim, td.rg, psd);
- ga_instruction_set::region_mim_instructions &rmi
- = gis.whole_instructions[rm];
+ auto &rmi = gis.all_instructions[rm];
rmi.m = td.m;
rmi.im = td.mim;
// rmi.interpolate_infos.clear();
@@ -6776,8 +6772,8 @@ namespace getfem {
// cout << "compilation finished "; ga_print_node(root, cout);
// cout << endl;
- if (version > 0) { // Assignment OR interpolation
- if (!td.varname_interpolation.empty()) { // assignment
+ if (phase != ga_workspace::ASSEMBLY) { // Assignment/interpolation
+ if (!td.varname_interpolation.empty()) {
auto *imd
= workspace.associated_im_data(td.varname_interpolation);
auto &V = const_cast<model_real_plain_vector &>
@@ -6963,8 +6959,8 @@ namespace getfem {
void ga_function_exec(ga_instruction_set &gis) {
- for (auto &&instr : gis.whole_instructions) {
- ga_instruction_list &gil = instr.second.instructions;
+ for (auto &&instr : gis.all_instructions) {
+ const auto &gil = instr.second.instructions;
for (size_type j = 0; j < gil.size(); ++j) j += gil[j]->exec();
}
}
@@ -6978,16 +6974,16 @@ namespace getfem {
for (const std::string &t : gis.transformations)
workspace.interpolate_transformation(t)->init(workspace);
- for (auto &&instr : gis.whole_instructions) {
+ for (auto &&instr : gis.all_instructions) {
const getfem::mesh_im &mim = *(instr.first.mim());
const mesh_region ®ion = *(instr.first.region());
const getfem::mesh &m = *(instr.second.m);
GMM_ASSERT1(&m == &(gic.linked_mesh()),
"Incompatibility of meshes in interpolation");
- ga_instruction_list &gilb = instr.second.begin_instructions;
- ga_instruction_list &gile = instr.second.elt_instructions;
- ga_instruction_list &gil = instr.second.instructions;
+ const auto &gilb = instr.second.begin_instructions;
+ const auto &gile = instr.second.elt_instructions;
+ const auto &gil = instr.second.instructions;
// iteration on elements (or faces of elements)
std::vector<size_type> ind;
@@ -7059,30 +7055,6 @@ namespace getfem {
gic.finalize();
}
- 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) {
- gis.ctx = ctx_x;
- gis.Normal = Normal;
- gmm::clear(workspace.assembled_tensor().as_vector());
- gis.nbpt = 1;
- gis.ipt = 0;
- gis.pai = 0;
-
- for (auto &&instr : gis.whole_instructions) {
- const getfem::mesh &m = *(instr.second.m);
- GMM_ASSERT1(&m == &interp_mesh,
- "Incompatibility of meshes in interpolation");
- ga_instruction_list &gilb = instr.second.begin_instructions;
- for (size_type j = 0; j < gilb.size(); ++j) j += gilb[j]->exec();
- ga_instruction_list &gile = instr.second.elt_instructions;
- for (size_type j = 0; j < gile.size(); ++j) j+=gile[j]->exec();
- ga_instruction_list &gil = instr.second.instructions;
- for (size_type j = 0; j < gil.size(); ++j) j += gil[j]->exec();
- }
- }
-
void ga_exec(ga_instruction_set &gis, ga_workspace &workspace) {
base_matrix G1, G2;
base_small_vector un;
@@ -7091,7 +7063,7 @@ namespace getfem {
for (const std::string &t : gis.transformations)
workspace.interpolate_transformation(t)->init(workspace);
- for (auto &instr : gis.whole_instructions) {
+ for (auto &instr : gis.all_instructions) {
const getfem::mesh_im &mim = *(instr.first.mim());
psecondary_domain psd = instr.first.psd();
const getfem::mesh &m = *(instr.second.m);
diff --git a/src/getfem_generic_assembly_interpolation.cc
b/src/getfem_generic_assembly_interpolation.cc
index 7702ee8..9fdddf7 100644
--- a/src/getfem_generic_assembly_interpolation.cc
+++ b/src/getfem_generic_assembly_interpolation.cc
@@ -636,8 +636,23 @@ namespace getfem {
bool compute_derivatives) const {
int ret_type = 0;
- ga_interpolation_single_point_exec(local_gis, local_workspace, ctx_x,
- Normal, m);
+ local_gis.ctx = ctx_x;
+ local_gis.Normal = Normal;
+ local_gis.nbpt = 1;
+ local_gis.ipt = 0;
+ local_gis.pai = 0;
+ gmm::clear(local_workspace.assembled_tensor().as_vector());
+
+ for (auto &&instr : local_gis.all_instructions) {
+ GMM_ASSERT1(instr.second.m == &m,
+ "Incompatibility of meshes in interpolation");
+ auto &gilb = instr.second.begin_instructions;
+ for (size_type j = 0; j < gilb.size(); ++j) j += gilb[j]->exec();
+ auto &gile = instr.second.elt_instructions;
+ for (size_type j = 0; j < gile.size(); ++j) j+=gile[j]->exec();
+ auto &gil = instr.second.instructions;
+ for (size_type j = 0; j < gil.size(); ++j) j += gil[j]->exec();
+ }
GMM_ASSERT1(local_workspace.assembled_tensor().size()==target_mesh.dim(),
"Wrong dimension of the transformation expression");
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index efbc6d3..4a7237b 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -45,14 +45,14 @@ namespace getfem {
void ga_workspace::add_fem_variable
(const std::string &name, const mesh_fem &mf,
const gmm::sub_interval &I, const model_real_plain_vector &VV) {
- variables[name] = var_description(true, true, &mf, I, &VV, 0, 1);
+ variables.emplace(name, var_description(true, &mf, 0, I, &VV, 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)));
+ variables.emplace(name, var_description(true, 0, 0, I, &VV,
+ dim_type(gmm::vect_size(VV))));
}
void ga_workspace::add_fem_constant
@@ -62,22 +62,22 @@ namespace getfem {
<< "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);
+ variables.emplace(name, var_description(false, &mf, 0,
+ gmm::sub_interval(), &VV, 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));
+ variables.emplace(name, var_description(false, 0, 0,
+ gmm::sub_interval(), &VV,
+ 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()));
+ variables.emplace(name, var_description
+ (false, 0, &imd, gmm::sub_interval(), &VV,
+ gmm::vect_size(VV)/(imd.nb_filtered_index() * imd.nb_tensor_elem())));
}
bool ga_workspace::variable_exists(const std::string &name) const {
@@ -390,7 +390,7 @@ namespace getfem {
const mesh_im &mim, const mesh_region &rg,
const std::string &expr,
size_type add_derivative_order,
- bool function_expr, size_type for_interpolation,
+ bool function_expr, operation_type op_type,
const std::string varname_interpolation) {
if (tree.root) {
// Eliminate the term if it corresponds to disabled variables
@@ -407,8 +407,8 @@ namespace getfem {
bool remain = true;
size_type order = 0, ind_tree = 0;
- if (for_interpolation)
- order = size_type(-1) - add_derivative_order;
+ if (op_type != ga_workspace::ASSEMBLY)
+ order = add_derivative_order;
else {
switch(tree.root->test_function_type) {
case 0: order = 0; break;
@@ -430,7 +430,7 @@ namespace getfem {
td.name_test2 == tree.root->name_test2 &&
td.interpolate_name_test2 == tree.root->interpolate_name_test2 &&
td.rg == &rg &&
- td.interpolation == for_interpolation &&
+ td.operation == op_type &&
td.varname_interpolation == varname_interpolation) {
ga_tree &ftree = *(td.ptree);
@@ -461,11 +461,11 @@ namespace getfem {
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().operation = op_type;
trees.back().varname_interpolation = varname_interpolation;
}
- if (for_interpolation == 0 && order < add_derivative_order) {
+ if (op_type == ga_workspace::ASSEMBLY && 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);
@@ -483,7 +483,7 @@ namespace getfem {
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);
+ function_expr, op_type, varname_interpolation);
}
}
}
@@ -539,7 +539,7 @@ namespace getfem {
// cout << "adding tree " << ga_tree_to_string(ltree) << endl;
max_order = std::max(ltree.root->nb_test_functions(), max_order);
add_tree(ltree, mim.linked_mesh(), mim, rg, expr,
- add_derivative_order, true, 0, "");
+ add_derivative_order, true);
}
}
}
@@ -555,7 +555,7 @@ namespace getfem {
// 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, "");
+ expr, 0, true);
}
}
@@ -570,7 +570,8 @@ namespace getfem {
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, "");
+ add_tree(tree, m, dummy_mesh_im(), rg, expr, 0, false,
+ ga_workspace::PRE_ASSIGNMENT);
}
}
@@ -586,7 +587,8 @@ namespace getfem {
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, "");
+ add_tree(tree, m, mim, rg, expr, 0, false,
+ ga_workspace::PRE_ASSIGNMENT);
}
}
@@ -604,7 +606,9 @@ namespace getfem {
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),
+ add_tree(tree, m, mim, rg, expr, order+1, false,
+ before ? ga_workspace::PRE_ASSIGNMENT
+ : ga_workspace::POST_ASSIGNMENT,
varname);
}
}
@@ -873,7 +877,7 @@ namespace getfem {
void ga_workspace::tree_description::copy(const tree_description& td) {
order = td.order;
- interpolation = td.interpolation;
+ operation = td.operation;
varname_interpolation = td.varname_interpolation;
name_test1 = td.name_test1;
name_test2 = td.name_test2;