[Getfem-commits] [getfem-commits] branch devel-logari81-internal-variables updated: Fix in proccessing of multiple coupled internal variables
This is an automated email from the git hooks/post-receive script.
logari81 pushed a commit to branch devel-logari81-internal-variables
in repository getfem.
The following commit(s) were added to
refs/heads/devel-logari81-internal-variables by this push:
new 19ccc74 Fix in proccessing of multiple coupled internal variables
19ccc74 is described below
commit 19ccc748b9dd0dbc50391074d1e3eeceb8f93d53
Author: Konstantinos Poulios
AuthorDate: Wed Jan 29 22:08:50 2020 +0100
Fix in proccessing of multiple coupled internal variables
---
src/getfem_generic_assembly_compile_and_exec.cc | 134 +---
1 file changed, 73 insertions(+), 61 deletions(-)
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index fe64993..f033698 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -5055,17 +5055,16 @@ namespace getfem {
std::vector RQloc;
base_tensor invK, Kqqjj;
base_vector Rqq;
-std::vector partQ, partJ;
-size_type Qsize, Jsize;
+std::vector> partQ, partJ;
virtual int exec() { // implementation can be optimized
GA_DEBUG_INFO("Instruction: variable cluster subdiagonal condensation");
// copy from KQQ to invK
- for (size_type q1=0; q1 < Qsize; ++q1) {
-size_type qq1start = partQ[q1], qq1end = partQ[q1+1];
-for (size_type q2=0; q2 < Qsize; ++q2) {
+ for (const auto : partQ) {
+size_type q1 = qqq1[0], qq1start = qqq1[1], qq1end = qqq1[2];
+for (const auto : partQ) {
+ size_type q2 = qqq2[0], qq2start = qqq2[1], qq2end = qqq2[2];
if (KQQloc(q1,q2)) {
auto itr = KQQloc(q1,q2)->cbegin();
-size_type qq2start = partQ[q2], qq2end = partQ[q2+1];
GMM_ASSERT1(KQQloc(q1,q2)->size()
== (qq1end-qq1start)*(qq2end-qq2start),
"Internal error");
@@ -5079,40 +5078,44 @@ namespace getfem {
bgeot::lu_inverse(&(invK[0]), invK.size(0));
// Resize Kqqjj as primary variable sizes may change dynamically
- partJ.clear();
- partJ.resize(Jsize,0);
- for (size_type j=0; j < Jsize; ++j)
-for (size_type q=0; q < Qsize; ++q)
+ size_type prev_j(0);
+ for (auto & : partJ) {
+size_type j=jjj[0];
+size_type new_j(0);
+for (const auto : partQ) {
+ size_type q=qqq[0];
if (KQJloc(q,j)) {
-if (partJ[j] == 0)
- partJ[j] = KQJloc(q,j)->size(1);
-else
- GMM_ASSERT1(partJ[j] == KQJloc(q,j)->size(1), "Internal error");
+if (new_j) {
+ GMM_ASSERT1(new_j == KQJloc(q,j)->size(1), "Internal error");
+} else
+ new_j = KQJloc(q,j)->size(1);
}
- for (size_type jj=1; jj < partJ.size(); ++jj)
-partJ[jj] += partJ[jj-1];
- partJ.insert(partJ.begin(), 0);
+}
+// Resize KQJprime submatrices to match KQJloc sizes
+for (const auto : partQ) {
+ size_type q=qqq[0];
+ KQJprime(q,j)->adjust_sizes(qqq[2]-qqq[1], new_j);
+}
+jjj[1] = prev_j;
+prev_j += new_j;
+jjj[2] = prev_j;
+ }
- Kqqjj.adjust_sizes(partQ.back(), partJ.back());
+ Kqqjj.adjust_sizes(partQ.back()[2], partJ.back()[2]);
gmm::clear(Kqqjj.as_vector());
gmm::clear(Rqq);
- // Resize KQJprime submatrices to match KQJloc sizes
- for (size_type j=0; j < Jsize; ++j)
-for (size_type q=0; q < Qsize; ++q)
- KQJprime(q,j)->adjust_sizes(partQ[q+1]-partQ[q],
partJ[j+1]-partJ[j]);
-
// multiply invK with all submatrices in KQJloc and RQloc and store
// the results in Kqqjj and Rqq
- for (size_type j=0; j < Jsize; ++j) {
-size_type jjstart = partJ[j], jjend = partJ[j+1];
-for (size_type q2=0; q2 < Qsize; ++q2) {
+ for (const auto : partJ) {
+size_type j = jjj[0], jjstart = jjj[1], jjend = jjj[2];
+for (const auto : partQ) {
+ size_type q2 = qqq2[0], qq2start = qqq2[1], qq2end = qqq2[2];
if (KQJloc(q2,j)) {
auto itr = KQJloc(q2,j)->begin(); // auto = KQJloc(q2,j);
-size_type qqstart = partQ[q2], qqend = partQ[q2+1];
for (size_type jj=jjstart; jj < jjend; ++jj) {
- for (size_type qq2=qqstart; qq2 < qqend; ++qq2, ++itr) {
-for (size_type qq1=0; qq1 < partQ.back(); ++qq1) {
+ for (size_type qq2=qq2start; qq2 < qq2end; ++qq2, ++itr) {
+for (size_type qq1=0; qq1 < partQ.back()[2]; ++qq1) {
Kqqjj(qq1,jj) += invK(qq1,qq2)*(*itr);
// Kqqjj(qq1,jj) +=
invKqq(qq1,qq2)*mat(qq2-qqstart,jj-jjstart);
} // for qq1
@@ -5120,30 +5123,31 @@ namespace getfem {
} // for jj
[Getfem-commits] (no subject)
branch: devel-logari81-internal-variables
commit 78b191f942d27d95da95d6dec284c0f051531cb1
Author: Konstantinos Poulios
AuthorDate: Wed Jan 29 16:47:37 2020 +0100
Minor changes
---
src/getfem/getfem_models.h | 16 +---
src/getfem_generic_assembly_compile_and_exec.cc | 12 +++-
src/getfem_models.cc| 16
3 files changed, 24 insertions(+), 20 deletions(-)
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 3e0657d..62ec9ee 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -271,12 +271,14 @@ namespace getfem {
typedef std::vector termlist;
-enum build_version { BUILD_RHS = 1,
- BUILD_MATRIX = 2,
- BUILD_ALL = 3,
- BUILD_ON_DATA_CHANGE = 4,
- BUILD_WITH_COMPLETE_RHS = 8,
- BUILD_COMPLETE_RHS = 9, };
+enum build_version {
+ BUILD_RHS = 1,
+ BUILD_MATRIX = 2,
+ BUILD_ALL = 3,
+ BUILD_ON_DATA_CHANGE = 4,
+ BUILD_WITH_LIN = 8, // forced calculation of linear terms
+ BUILD_RHS_WITH_LIN = 9, // = BUILD_RHS | BUILD_WITH_LIN_RHS
+};
protected:
@@ -356,7 +358,7 @@ namespace getfem {
std::string secondary_domain;
gen_expr(const std::string _, const mesh_im _,
size_type region_, const std::string )
- : expr(expr_), mim(mim_), region(region_), secondary_domain(secdom) {}
+: expr(expr_), mim(mim_), region(region_), secondary_domain(secdom) {}
};
// Structure for assignment in assembly
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 11c2a72..fe64993 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -5445,9 +5445,11 @@ namespace getfem {
}
cout << "sub_tree_are_equal = " << int(sub_tree_are_equal(pnode,
pnode1, workspace, 1)) << endl;
std::stringstream ss;
-ss << "Detected wrong equivalent nodes: \n";
-ga_print_node(pnode, ss); ss << "\n and \n"; ga_print_node(pnode1, ss);
-ss << "\nNo problem, but hash values could be adapted." << endl;
+ss << "Detected wrong equivalent nodes:" << endl;
+ga_print_node(pnode, ss);
+ss << endl << " and " << endl;
+ga_print_node(pnode1, ss);
+ss << endl << "No problem, but hash values could be adapted." << endl;
GMM_TRACE2(ss.str());
}
}
@@ -8019,13 +8021,13 @@ namespace getfem {
(Kq1j2pr, KQJpr, gis.ctx, gis.ctx,
I1, imd1, alpha1, I2, imd2, alpha2, gis.ipt); //
constructor without gis.coeff
rmi.instructions.push_back(std::move(pgai));
- }
+ } // for j2
const bool initialize = true;
pgai = std::make_shared
(*(CC.RQpr[q1]), workspace.assembled_vector(), // <-
overwriting internal variables residual with internal solution
gis.ctx, I1, *imd1, gis.coeff, gis.ipt, initialize);
rmi.instructions.push_back(std::move(pgai));
-}
+} // for q1
}
// Add superdiagonal condensation instructions
diff --git a/src/getfem_models.cc b/src/getfem_models.cc
index 25e67c1..c8a5f1a 100644
--- a/src/getfem_models.cc
+++ b/src/getfem_models.cc
@@ -2391,7 +2391,7 @@ namespace getfem {
}
}
if (term.is_matrix_term && brick.pbr->is_linear()
-&& (!is_linear() || (version & BUILD_WITH_COMPLETE_RHS))) {
+&& (!is_linear() || (version & BUILD_WITH_LIN))) {
if (!(var1->is_disabled))
gmm::mult_add(brick.cmatlist[j],
gmm::scaled(var2->complex_value[0],
@@ -2411,7 +2411,7 @@ namespace getfem {
gmm::sub_vector(crhs, I2));
}
if (brick.pbr->is_linear()
- && (!is_linear() || (version & BUILD_WITH_COMPLETE_RHS))) {
+ && (!is_linear() || (version & BUILD_WITH_LIN))) {
gmm::mult_add(gmm::conjugated(brick.cmatlist[j]),
gmm::scaled(var1->complex_value[0],
complex_type(-alpha2)),
@@ -2459,7 +2459,7 @@ namespace getfem {
}
}
if (term.is_matrix_term && brick.pbr->is_linear()
-&& (!is_linear() || (version & BUILD_WITH_COMPLETE_RHS))) {
+&& (!is_linear() || (version & BUILD_WITH_LIN))) {
if (!(var1->is_disabled))
gmm::mult_add(brick.rmatlist[j],
gmm::scaled(var2->complex_value[0],
@@ -2479,7 +2479,7 @@ namespace getfem {
[Getfem-commits] (no subject)
branch: devel-logari81-internal-variables
commit 50ebaab6f372cc9435066732a13a9cb14498caa5
Author: Konstantinos Poulios
AuthorDate: Wed Jan 29 16:49:18 2020 +0100
Add support for internal variable condensation to the model object
---
src/getfem/getfem_accumulated_distro.h | 6 +-
src/getfem/getfem_models.h | 42 ++--
src/getfem_generic_assembly_workspace.cc | 4 +-
src/getfem_model_solvers.cc | 90 +++-
src/getfem_models.cc | 178 +++
5 files changed, 261 insertions(+), 59 deletions(-)
diff --git a/src/getfem/getfem_accumulated_distro.h
b/src/getfem/getfem_accumulated_distro.h
index ebf8226..a428202 100644
--- a/src/getfem/getfem_accumulated_distro.h
+++ b/src/getfem/getfem_accumulated_distro.h
@@ -173,12 +173,16 @@ namespace detail {
}
}
-operator T&(){
+T& get(){
if (distributed.num_threads() == 1 ||
distributed.this_thread() == 0) return original;
else return distributed;
}
+operator T&(){ // implicit conversion
+ return get();
+}
+
T& operator = (const T ){
return distributed = x;
}
diff --git a/src/getfem/getfem_models.h b/src/getfem/getfem_models.h
index 62ec9ee..813bb69 100644
--- a/src/getfem/getfem_models.h
+++ b/src/getfem/getfem_models.h
@@ -121,9 +121,14 @@ namespace getfem {
bool is_linear_;
bool is_symmetric_;
bool is_coercive_;
-mutable model_real_sparse_matrix rTM;// tangent matrix, real version
+mutable model_real_sparse_matrix
+ rTM, // tangent matrix (only primary variables), real version
+ internal_rTM; // coupling matrix between internal and primary vars (no
empty rows)
mutable model_complex_sparse_matrix cTM; // tangent matrix, complex version
-mutable model_real_plain_vector rrhs;
+mutable model_real_plain_vector
+ rrhs, // residual vector of primary variables (after
condensation)
+ full_rrhs,// residual vector of primary and internal variables
(pre-condensation)
+ internal_sol; // partial solution for internal variables (after
condensation)
mutable model_complex_plain_vector crhs;
mutable bool act_size_to_be_done;
dim_type leading_dim;
@@ -278,6 +283,10 @@ namespace getfem {
BUILD_ON_DATA_CHANGE = 4,
BUILD_WITH_LIN = 8, // forced calculation of linear terms
BUILD_RHS_WITH_LIN = 9, // = BUILD_RHS | BUILD_WITH_LIN_RHS
+ BUILD_WITH_INTERNAL = 16,
+ BUILD_RHS_WITH_INTERNAL = 17, // = BUILD_RHS | BUILD_WITH_INTERNAL
+ BUILD_MATRIX_CONDENSED = 18, // = BUILD_MATRIX | BUILD_WITH_INTERNAL
+ BUILD_ALL_CONDENSED = 19, // = BUILD_ALL | BUILD_WITH_INTERNAL
};
protected:
@@ -571,7 +580,13 @@ namespace getfem {
bool is_linear() const { return is_linear_; }
/** Total number of degrees of freedom in the model. */
-size_type nb_dof() const;
+size_type nb_dof(bool with_internal=false) const;
+
+/** Number of internal degrees of freedom in the model. */
+size_type nb_internal_dof() const { return nb_dof(true)-nb_dof(); }
+
+/** Number of primary degrees of freedom in the model. */
+size_type nb_primary_dof() const { return nb_dof(); }
/** Leading dimension of the meshes used in the model. */
dim_type leading_dimension() const { return leading_dim; }
@@ -898,10 +913,11 @@ namespace getfem {
size_type qdim_of_variable(const std::string ) const;
/** Gives the access to the tangent matrix. For the real version. */
-const model_real_sparse_matrix _tangent_matrix() const {
+const model_real_sparse_matrix &
+real_tangent_matrix(bool internal=false) const {
GMM_ASSERT1(!complex_version, "This model is a complex one");
context_check(); if (act_size_to_be_done) actualize_sizes();
- return rTM;
+ return internal ? internal_rTM : rTM;
}
/** Gives the access to the tangent matrix. For the complex version. */
@@ -913,19 +929,27 @@ namespace getfem {
/** Gives access to the right hand side of the tangent linear system.
For the real version. An assembly of the rhs has to be done first. */
-const model_real_plain_vector _rhs() const {
+const model_real_plain_vector _rhs(bool with_internal=false) const {
GMM_ASSERT1(!complex_version, "This model is a complex one");
context_check(); if (act_size_to_be_done) actualize_sizes();
- return rrhs;
+ return (with_internal && gmm::vect_size(full_rrhs)) ? full_rrhs : rrhs;
}
/** Gives write access to the right hand side of the tangent linear system.
Some solvers need to manipulate the model rhs directly so that for
example internal condensed variables can be treated properly. */
-model_real_plain_vector _real_rhs() const {
+model_real_plain_vector _real_rhs(bool with_internal=false) const {
+ GMM_ASSERT1(!complex_version,
[Getfem-commits] (no subject)
branch: devel-logari81-internal-variables
commit 232385db09a7ff588cb8144fd1e8aa79ce6a32b6
Author: Konstantinos Poulios
AuthorDate: Wed Jan 29 09:47:03 2020 +0100
Fix size of internal variable condensation matrix
- Including empty rows for keeping the implementation simpler
---
src/getfem/getfem_generic_assembly.h| 4 +++-
src/getfem_generic_assembly_compile_and_exec.cc | 7 +++
src/getfem_generic_assembly_workspace.cc| 6 ++
3 files changed, 12 insertions(+), 5 deletions(-)
diff --git a/src/getfem/getfem_generic_assembly.h
b/src/getfem/getfem_generic_assembly.h
index 16b77b6..367d227 100644
--- a/src/getfem/getfem_generic_assembly.h
+++ b/src/getfem/getfem_generic_assembly.h
@@ -413,7 +413,9 @@ namespace getfem {
KQJpr = std::shared_ptr
(std::shared_ptr(), _); // alias
}
-// getter functions for condensation matrix/vectors
+/** getter function for condensation matrix
+Its size is (nb_primary_dof()+nb_internal_dof()) x nb_primary_dof()
+but its first nb_primary_dof() rows are empty */
const model_real_sparse_matrix _coupling_matrix() const
{ return *KQJpr; }
model_real_sparse_matrix _coupling_matrix() { return *KQJpr; }
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 5c301ba..11c2a72 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -8000,6 +8000,13 @@ namespace getfem {
const base_tensor = *(CC.KQJpr(q1,j2)); // <- input
model_real_sparse_matrix
= workspace.internal_coupling_matrix(); // <- output
+GMM_ASSERT1(
+ gmm::mat_nrows(KQJpr) == workspace.nb_primary_dof()
+ +workspace.nb_internal_dof() &&
+ gmm::mat_ncols(KQJpr) == workspace.nb_primary_dof(),
+ "The internal coupling matrix needs to be allocated with "
+ "nb_primary_dof()+nb_internal_dof() number of rows, even if "
+ "only the last nb_internal_dof() rows are filled in.");
if (mf2)
pgai =
std::make_shared
diff --git a/src/getfem_generic_assembly_workspace.cc
b/src/getfem_generic_assembly_workspace.cc
index 50c0db1..4b7517f 100644
--- a/src/getfem_generic_assembly_workspace.cc
+++ b/src/getfem_generic_assembly_workspace.cc
@@ -789,9 +789,9 @@ namespace getfem {
// gmm::mat_ncols(*K) == nb_prim_dof, "Wrong sizes");
if (KQJpr.use_count()) {
gmm::clear(*KQJpr);
-gmm::resize(*KQJpr, nb_intern_dof, nb_prim_dof); // redundant if
condensation == false
+gmm::resize(*KQJpr, nb_prim_dof+nb_intern_dof, nb_prim_dof); //
redundant if condensation == false
} else if (condensation)
-GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_intern_dof &&
+GMM_ASSERT1(gmm::mat_nrows(*KQJpr) == nb_prim_dof+nb_intern_dof &&
gmm::mat_ncols(*KQJpr) == nb_prim_dof, "Wrong sizes");
gmm::clear(col_unreduced_K);
gmm::clear(row_unreduced_K);
@@ -894,8 +894,6 @@ namespace getfem {
GMM_ASSERT1(I1.last() <= nb_prim_dof, "Internal error");
gmm::add(M, gmm::sub_matrix(*K, I1, I2)); // -> *K
} else { // vname1 is an internal variable
- I1.min -= first_internal_dof();
- I1.max -= first_internal_dof();
gmm::add(M, gmm::sub_matrix(*KQJpr, I1, I2)); // ->
*KQJpr
}
}
[Getfem-commits] (no subject)
branch: devel-logari81-internal-variables
commit ff5eb9c9271cb82fd96b6044728bae6173203c81
Author: Konstantinos Poulios
AuthorDate: Wed Jan 29 08:57:05 2020 +0100
Fix wrong scaling of coupling matrix in internal variable condensation
---
src/getfem_generic_assembly_compile_and_exec.cc | 29 +++--
1 file changed, 27 insertions(+), 2 deletions(-)
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 890b1b5..5c301ba 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -4534,6 +4534,7 @@ namespace getfem {
protected:
const bool false_=false;
const size_type zero_=0;
+const scalar_type one_=1;
};
@@ -4714,6 +4715,18 @@ namespace getfem {
(t_, ctx1_, ctx2_, a1, a2, coeff_, zero_, ipt_, false),
Kxr(K_), Kxu(K_), I1(_), I2__(_), I2(I2__),
imd1(imd1_), mf2__(_), mf2(mf2__), reduced_mf2(false_) {}
+
+ga_instruction_matrix_assembly_imd_mf // constructor without coeff (== 1)
+(const base_tensor _, model_real_sparse_matrix _,
+ const fem_interpolation_context _,
+ const fem_interpolation_context _,
+ const gmm::sub_interval _, const im_data *imd1_, const scalar_type ,
+ const gmm::sub_interval _, const mesh_fem _, const scalar_type ,
+ const size_type _)
+ : ga_instruction_matrix_assembly_base
+(t_, ctx1_, ctx2_, a1, a2, one_, zero_, ipt_, false),
+Kxr(K_), Kxu(K_), I1(_), I2__(_), I2(I2__),
+imd1(imd1_), mf2__(_), mf2(mf2__), reduced_mf2(false_) {}
};
struct ga_instruction_matrix_assembly_mf_imd
@@ -4816,6 +4829,17 @@ namespace getfem {
: ga_instruction_matrix_assembly_base
(t_, ctx1_, ctx2_, a1, a2, coeff_, zero_, ipt_, false),
K(K_), I1(I1_), I2(I2_), imd1(imd1_), imd2(imd2_) {}
+
+ga_instruction_matrix_assembly_imd_imd // constructor without coeff (== 1)
+(const base_tensor _, model_real_sparse_matrix _,
+ const fem_interpolation_context _,
+ const fem_interpolation_context _,
+ const gmm::sub_interval _, const im_data *imd1_, const scalar_type ,
+ const gmm::sub_interval _, const im_data *imd2_, const scalar_type ,
+ const size_type _)
+ : ga_instruction_matrix_assembly_base
+(t_, ctx1_, ctx2_, a1, a2, one_, zero_, ipt_, false),
+K(K_), I1(I1_), I2(I2_), imd1(imd1_), imd2(imd2_) {}
};
@@ -7980,12 +8004,13 @@ namespace getfem {
pgai =
std::make_shared
(Kq1j2pr, KQJpr, gis.ctx, gis.ctx,
- I1, imd1, alpha1, I2, *mf2, alpha2, gis.coeff, gis.ipt);
// TODO: name_test2 variable group
+ I1, imd1, alpha1, I2, *mf2, alpha2, gis.ipt); //
constructor without gis.coeff
+// TODO: name_test2 variable group
else // for global variable imd2 == 0
pgai =
std::make_shared
(Kq1j2pr, KQJpr, gis.ctx, gis.ctx,
- I1, imd1, alpha1, I2, imd2, alpha2, gis.coeff, gis.ipt);
+ I1, imd1, alpha1, I2, imd2, alpha2, gis.ipt); //
constructor without gis.coeff
rmi.instructions.push_back(std::move(pgai));
}
const bool initialize = true;
