[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: |
Thu, 20 Dec 2018 03:46:29 -0500 (EST) |
branch: integration-point-variables
commit b300253d05c063ab1125beff4d948196b2b80af1
Author: Konstantinos Poulios <address@hidden>
Date: Thu Dec 20 09:46:12 2018 +0100
Enable integration point (internal) variables
---
src/getfem/getfem_generic_assembly.h | 13 ++--
src/getfem/getfem_models.h | 4 ++
src/getfem_generic_assembly_compile_and_exec.cc | 72 +++++++++++++++++---
src/getfem_generic_assembly_semantic.cc | 87 ++++++++++++++++++-------
src/getfem_generic_assembly_workspace.cc | 6 ++
src/getfem_models.cc | 11 ++++
6 files changed, 154 insertions(+), 39 deletions(-)
diff --git a/src/getfem/getfem_generic_assembly.h
b/src/getfem/getfem_generic_assembly.h
index 9468406..fad863b 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -276,7 +276,8 @@ namespace getfem {
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,10 +287,9 @@ namespace getfem {
return q;
}
- var_description(bool is_var, bool is_fem,
- const mesh_fem *mf_, gmm::sub_interval I_,
- const model_real_plain_vector *v, const im_data *imd_,
- size_type Q)
+ var_description(bool is_var, bool is_fem, const mesh_fem *mf_,
+ gmm::sub_interval I_, const model_real_plain_vector *v,
+ const im_data *imd_, size_type Q)
: is_variable(is_var), is_fem_dofs(is_fem), mf(mf_), I(I_), V(v),
imd(imd_), qdims(1)
{
@@ -439,6 +439,9 @@ namespace getfem {
void add_fem_variable(const std::string &name, const mesh_fem &mf,
const gmm::sub_interval &I,
const model_real_plain_vector &VV);
+ void add_im_variable(const std::string &name, const im_data &imd,
+ const gmm::sub_interval &I,
+ const model_real_plain_vector &VV);
void add_fixed_size_variable(const std::string &name,
const gmm::sub_interval &I,
const model_real_plain_vector &VV);
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 1b71bcc..82b35c9 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -750,6 +750,10 @@ namespace getfem {
}
+ /** Add (internal) variable, defined at integration points.*/
+ void add_im_variable(const std::string &name, const im_data &im_data,
+ size_type niter = 1);
+
/** Add data, defined at integration points.*/
void add_im_data(const std::string &name, const im_data &im_data,
size_type niter = 1);
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 243d409..abc7e9e 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -4075,6 +4075,32 @@ namespace getfem {
interpolate(interpolate_) {}
};
+ struct ga_instruction_imd_vector_assembly : public ga_instruction {
+ const base_tensor &t;
+ base_vector &V;
+ const fem_interpolation_context &ctx;
+ const gmm::sub_interval &I;
+ const im_data *imd;
+ scalar_type &coeff;
+ const size_type &ipt;
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: vector term assembly for im_data variable");
+ size_type ifirst = I.first();
+ size_type i = t.size()
+ * imd->filtered_index_of_point(ctx.convex_num(), ipt);
+ GMM_ASSERT1(i+t.size() <= I.size(), "Internal error
"<<i<<"+"<<t.size()<<" <= "<<I.size());
+ for (const auto &val : t.as_vector())
+ V[ifirst+(i++)] += coeff*val;
+ return 0;
+ }
+ ga_instruction_imd_vector_assembly
+ (const base_tensor &t_, base_vector &V_,
+ const fem_interpolation_context &ctx_, const gmm::sub_interval &I_,
+ const im_data *imd_, scalar_type &coeff_, const size_type &ipt_)
+ : t(t_), V(V_), ctx(ctx_), I(I_), imd(imd_), coeff(coeff_), ipt(ipt_)
+ {}
+ };
+
struct ga_instruction_vector_assembly : public ga_instruction {
base_tensor &t;
base_vector &V;
@@ -4209,6 +4235,7 @@ namespace getfem {
const gmm::sub_interval &In1, &In2;
const mesh_fem *mfn1, *mfn2;
const mesh_fem **mfg1, **mfg2;
+ const im_data *imd1, *imd2;
const scalar_type &coeff, &alpha1, &alpha2;
const size_type &nbpt, &ipt;
base_vector elem;
@@ -4217,7 +4244,10 @@ namespace getfem {
virtual int exec() {
GA_DEBUG_INFO("Instruction: matrix term assembly");
bool empty_weight = (coeff == scalar_type(0));
- if (ipt == 0 || interpolate) {
+ bool initialize = (ipt == 0) || interpolate || imd1 || imd2;
+ bool finalize = (ipt == nbpt-1) || interpolate || imd1 || imd2;
+
+ if (initialize) {
if (empty_weight) elem.resize(0);
elem.resize(t.size());
if (!empty_weight) {
@@ -4242,7 +4272,8 @@ namespace getfem {
for (; it != ite;) *it++ += (*itt++) * e;
// gmm::add(gmm::scaled(t.as_vector(), coeff*alpha1*alpha2), elem);
}
- if (ipt == nbpt-1 || interpolate) {
+
+ if (finalize) {
const mesh_fem *pmf1 = mfg1 ? *mfg1 : mfn1;
const mesh_fem *pmf2 = mfg2 ? *mfg2 : mfn2;
bool reduced = (pmf1 && pmf1->is_reduced())
@@ -4260,7 +4291,13 @@ namespace getfem {
size_type cv2 = pmf2 ? ctx2.convex_num() : s2;
size_type N = 1;
- dofs1.assign(s1, I1.first());
+ size_type ifirst1 = I1.first(), ifirst2 = I2.first();
+ if (imd1)
+ ifirst1 += s1*imd1->filtered_index_of_point(ctx1.convex_num(), ipt);
+ if (imd2)
+ ifirst2 += s2*imd2->filtered_index_of_point(ctx2.convex_num(), ipt);
+
+ dofs1.assign(s1, ifirst1);
if (pmf1) {
if (!(ctx1.is_convex_num_valid())) return 0;
N = ctx1.N();
@@ -4280,16 +4317,16 @@ namespace getfem {
for (size_type i=0; i < s1; ++i) dofs1[i] += i;
if (pmf1 == pmf2 && cv1 == cv2) {
- if (I1.first() == I2.first()) {
+ if (ifirst1 == ifirst2) {
add_elem_matrix_(K, dofs1, dofs1, dofs1_sort, elem, ninf*1E-14, N);
} else {
dofs2.resize(dofs1.size());
for (size_type i = 0; i < dofs1.size(); ++i)
- dofs2[i] = dofs1[i] + I2.first() - I1.first();
+ dofs2[i] = dofs1[i] + ifirst2 - ifirst1;
add_elem_matrix_(K, dofs1, dofs2, dofs1_sort, elem, ninf*1E-14, N);
}
} else {
- dofs2.assign(s2, I2.first());
+ dofs2.assign(s2, ifirst2);
if (pmf2) {
if (!(ctx2.is_convex_num_valid())) return 0;
N = std::max(N, ctx2.N());
@@ -4319,14 +4356,15 @@ namespace getfem {
const fem_interpolation_context &ctx2_,
const gmm::sub_interval &Ir1_, const gmm::sub_interval &In1_,
const gmm::sub_interval &Ir2_, const gmm::sub_interval &In2_,
- const mesh_fem *mfn1_, const mesh_fem **mfg1_,
- const mesh_fem *mfn2_, const mesh_fem **mfg2_,
+ const mesh_fem *mfn1_, const mesh_fem **mfg1_, const im_data *imd1_,
+ const mesh_fem *mfn2_, const mesh_fem **mfg2_, const im_data *imd2_,
const scalar_type &coeff_,
const scalar_type &alpha2_, const scalar_type &alpha1_,
const size_type &nbpt_, const size_type &ipt_, bool interpolate_)
: t(t_), Kr(Kr_), Kn(Kn_), ctx1(ctx1_), ctx2(ctx2_),
Ir1(Ir1_), Ir2(Ir2_), In1(In1_), In2(In2_),
mfn1(mfn1_), mfn2(mfn2_), mfg1(mfg1_), mfg2(mfg2_),
+ imd1(imd1_), imd2(imd2_),
coeff(coeff_), alpha1(alpha1_), alpha2(alpha2_),
nbpt(nbpt_), ipt(ipt_), interpolate(interpolate_),
dofs1(0), dofs2(0) {}
@@ -6761,7 +6799,10 @@ namespace getfem {
GMM_ASSERT1(root->tensor_proper_size() == 1,
"Invalid vector or tensor quantity. An order 1 "
"weak form has to be a scalar quantity");
- const mesh_fem *mf=workspace.associated_mf(root->name_test1);
+ const mesh_fem *mf
+ = workspace.associated_mf(root->name_test1);
+ const im_data *imd
+ = workspace.associated_im_data(root->name_test1);
add_interval_to_gis(workspace, root->name_test1, gis);
if (mf) {
@@ -6795,6 +6836,14 @@ namespace getfem {
(root->tensor(), workspace.unreduced_vector(),
workspace.assembled_vector(), ctx, *Ir, *In, mf, mfg,
gis.coeff, gis.nbpt, gis.ipt, interpolate);
+ } else if (imd) {
+ GMM_ASSERT1(root->interpolate_name_test1.size() == 0,
+ "Interpolate transformation on integration "
+ "point variable");
+ pgai = std::make_shared<ga_instruction_imd_vector_assembly>
+ (root->tensor(), workspace.assembled_vector(), gis.ctx,
+ workspace.interval_of_variable(root->name_test1),
+ imd, gis.coeff, gis.ipt);
} else {
pgai = std::make_shared<ga_instruction_vector_assembly>
(root->tensor(), workspace.assembled_vector(),
@@ -6811,6 +6860,9 @@ namespace getfem {
const mesh_fem
*mf1=workspace.associated_mf(root->name_test1);
const mesh_fem
*mf2=workspace.associated_mf(root->name_test2);
const mesh_fem **mfg1 = 0, **mfg2 = 0;
+ const im_data
+ *imd1 = workspace.associated_im_data(root->name_test1),
+ *imd2 = workspace.associated_im_data(root->name_test2);
const std::string &intn1 = root->interpolate_name_test1;
const std::string &intn2 = root->interpolate_name_test2;
bool secondary1 = intn1.size() &&
@@ -6902,7 +6954,7 @@ namespace getfem {
pgai = std::make_shared<ga_instruction_matrix_assembly<>>
(root->tensor(), workspace.unreduced_matrix(),
workspace.assembled_matrix(), ctx1, ctx2,
- *Ir1, *In1, *Ir2, *In2, mf1, mfg1, mf2, mfg2,
+ *Ir1, *In1, *Ir2, *In2, mf1, mfg1, imd1, mf2, mfg2,
imd2,
gis.coeff, *alpha1, *alpha2, gis.nbpt, gis.ipt,
interpolate);
}
diff --git a/src/getfem_generic_assembly_semantic.cc
b/src/getfem_generic_assembly_semantic.cc
index baa422f..041235d 100644
--- a/src/getfem_generic_assembly_semantic.cc
+++ b/src/getfem_generic_assembly_semantic.cc
@@ -471,13 +471,15 @@ namespace getfem {
case GA_NODE_XFEM_MINUS_HESS_TEST: case GA_NODE_XFEM_MINUS_DIVERG_TEST:
{
const mesh_fem *mf = workspace.associated_mf(pnode->name);
+ const im_data *imd = workspace.associated_im_data(pnode->name);
size_type t_type = pnode->test_function_type;
if (t_type == 1) {
pnode->name_test1 = pnode->name;
pnode->interpolate_name_test1 = pnode->interpolate_name;
pnode->interpolate_name_test2 = pnode->name_test2 = "";
- pnode->qdim1 = (mf ? workspace.qdim(pnode->name)
- : gmm::vect_size(workspace.value(pnode->name)));
+ pnode->qdim1 = (mf || imd)
+ ? workspace.qdim(pnode->name)
+ : gmm::vect_size(workspace.value(pnode->name));
if (option == 1)
workspace.test1.insert
(var_trans_pair(pnode->name_test1,
@@ -489,8 +491,9 @@ namespace getfem {
pnode->interpolate_name_test1 = pnode->name_test1 = "";
pnode->name_test2 = pnode->name;
pnode->interpolate_name_test2 = pnode->interpolate_name;
- pnode->qdim2 = (mf ? workspace.qdim(pnode->name)
- : gmm::vect_size(workspace.value(pnode->name)));
+ pnode->qdim2 = (mf || imd)
+ ? workspace.qdim(pnode->name)
+ : gmm::vect_size(workspace.value(pnode->name));
if (option == 1)
workspace.test2.insert
(var_trans_pair(pnode->name_test2,
@@ -499,11 +502,11 @@ namespace getfem {
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
}
- if (!mf) {
- size_type n = workspace.qdim(pnode->name);
- if (!n)
- ga_throw_error(pnode->expr, pnode->pos,
- "Invalid null size of variable");
+ size_type n = workspace.qdim(pnode->name);
+ if (!n)
+ ga_throw_error(pnode->expr, pnode->pos,
+ "Invalid null size of variable");
+ if (!mf & !imd) { // global variable
if (n == 1) {
pnode->init_vector_tensor(1);
pnode->tensor()[0] = scalar_type(1);
@@ -513,8 +516,22 @@ namespace getfem {
pnode->test_function_type = t_type;
for (size_type i = 0; i < n; ++i)
for (size_type j = 0; j < n; ++j)
- pnode->tensor()(i,j) = (i==j) ? scalar_type(1) :
scalar_type(0);
+ pnode->tensor()(i,j) = (i==j) ? scalar_type(1)
+ : scalar_type(0);
}
+ } else if (imd) {
+ bgeot::multi_index mii = workspace.qdims(pnode->name);
+ if (n == 1 && mii.size() <= 1) {
+ mii.resize(1);
+ mii[0] = 1;
+ } else
+ mii.insert(mii.begin(), n);
+ pnode->t.adjust_sizes(mii);
+ GMM_ASSERT1(pnode->tensor().size() == n*n, "Internal error");
+ auto itt = pnode->tensor().begin(); // set t equal to identity
+ for (size_type i = 0; i < n; ++i)
+ for (size_type j = 0; j < n; ++j)
+ *itt++ = (i == j) ? scalar_type(1) : scalar_type(0);
}
}
break;
@@ -1730,8 +1747,8 @@ namespace getfem {
workspace.test1.insert
(var_trans_pair(pnode->name_test1,
pnode->interpolate_name_test1));
- pnode->qdim1 = mf ? workspace.qdim(name)
- : gmm::vect_size(workspace.value(name));
+ pnode->qdim1 = (mf || imd) ? workspace.qdim(name)
+ : gmm::vect_size(workspace.value(name));
if (!(pnode->qdim1))
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
@@ -1742,14 +1759,14 @@ namespace getfem {
workspace.test2.insert
(var_trans_pair(pnode->name_test2,
pnode->interpolate_name_test2));
- pnode->qdim2 = mf ? workspace.qdim(name)
- : gmm::vect_size(workspace.value(name));
+ pnode->qdim2 = (mf || imd) ? workspace.qdim(name)
+ : gmm::vect_size(workspace.value(name));
if (!(pnode->qdim2))
ga_throw_error(pnode->expr, pnode->pos,
"Invalid null size of variable");
}
- if (!mf && (test || !imd)) {
+ if (!mf && !imd) { // global variable
if (prefix_id)
ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
"Divergence cannot be evaluated for fixed size data.");
@@ -1773,20 +1790,44 @@ namespace getfem {
pnode->init_matrix_tensor(n,n);
for (size_type i = 0; i < n; ++i)
for (size_type j = 0; j < n; ++j)
- pnode->tensor()(i,j)
- = ((i == j) ? scalar_type(1) : scalar_type(0));
+ pnode->tensor()(i,j) = (i == j) ? scalar_type(1)
+ : scalar_type(0);
} else {
pnode->t.adjust_sizes(workspace.qdims(name));
gmm::copy(workspace.value(name), pnode->tensor().as_vector());
}
}
- } else if (!test && imd) {
+ } else if (imd) { // im_data variable
+ size_type q = workspace.qdim(name);
+ bgeot::multi_index mii = workspace.qdims(name);
+
+ if (!q) ga_throw_error(pnode->expr, pnode->pos,
+ "Invalid null size of variable " << name);
+ if (mii.size() > 6)
+ ga_throw_error(pnode->expr, pnode->pos,
+ "Tensor with too many dimensions. Limited to 6");
if (prefix_id)
ga_throw_error(pnode->expr, pnode->pos, "Gradient, Hessian or "
"Divergence cannot be evaluated for im data.");
- pnode->node_type = GA_NODE_VAL;
- pnode->t.adjust_sizes(workspace.qdims(name));
- } else {
+
+ pnode->node_type = test ? GA_NODE_VAL_TEST : GA_NODE_VAL;
+
+ if (test) {
+ if (q == 1 && mii.size() <= 1) {
+ mii.resize(1);
+ mii[0] = 1;
+ } else
+ mii.insert(mii.begin(), q);
+ }
+ pnode->t.adjust_sizes(mii);
+ if (test) {
+ GMM_ASSERT1(pnode->tensor().size() == q*q, "Internal error");
+ auto itt = pnode->tensor().begin(); // set t equal to identity
+ for (size_type i = 0; i < q; ++i)
+ for (size_type j = 0; j < q; ++j)
+ *itt++ = (i == j) ? scalar_type(1) : scalar_type(0);
+ }
+ } else { // mesh_fem variable
size_type q = workspace.qdim(name);
size_type n = mf->linked_mesh().dim();
bgeot::multi_index mii = workspace.qdims(name);
@@ -1805,10 +1846,8 @@ namespace getfem {
if (test && q == 1 && mii.size() <= 1) {
mii.resize(1);
mii[0] = 2;
- } else if (test) {
+ } else if (test)
mii.insert(mii.begin(), 2);
- pnode->t.adjust_sizes(mii);
- }
break;
case 1: // grad
pnode->node_type = test ? GA_NODE_GRAD_TEST : GA_NODE_GRAD;
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index 77f6a5c..a06aeef 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -48,6 +48,12 @@ namespace getfem {
variables[name] = var_description(true, true, &mf, I, &VV, 0, 1);
}
+ void ga_workspace::add_im_variable
+ (const std::string &name, const im_data &imd,
+ const gmm::sub_interval &I, const model_real_plain_vector &VV) {
+ variables[name] = var_description(true, false, 0, I, &VV, &imd, 1);
+ }
+
void ga_workspace::add_fixed_size_variable
(const std::string &name,
const gmm::sub_interval &I, const model_real_plain_vector &VV) {
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index 1f7c6e1..ad29100 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -842,6 +842,17 @@ namespace getfem {
gmm::copy(t.as_vector(), this->set_complex_variable(name));
}
+ void model::add_im_variable(const std::string &name,
+ const im_data &im_data,
+ size_type niter) {
+ check_name_validity(name);
+ variables[name] = var_description(true, is_complex(), false, niter);
+ variables[name].pim_data = &im_data;
+ variables[name].set_size();
+ add_dependency(im_data);
+ act_size_to_be_done = true;
+ }
+
void model::add_im_data(const std::string &name, const im_data &im_data,
size_type niter) {
check_name_validity(name);