[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: |
Fri, 24 Aug 2018 04:25:20 -0400 (EDT) |
branch: master
commit 05a4ea6737507abd5e695248bd677813092f10c1
Author: Yves Renard <address@hidden>
Date: Fri Aug 24 10:24:56 2018 +0200
A mean to select the output matrix/vector part corresponding to a variable
in the generica ssembly of gf_asm
---
interface/src/gf_asm.cc | 71 ++++++++++++++++------
.../tests/matlab/demo_Nitsche_contact_by_hand.m | 6 +-
src/getfem_generic_assembly_compile_and_exec.cc | 27 ++++++--
3 files changed, 79 insertions(+), 25 deletions(-)
diff --git a/interface/src/gf_asm.cc b/interface/src/gf_asm.cc
index 5792ecf..1e76be1 100644
--- a/interface/src/gf_asm.cc
+++ b/interface/src/gf_asm.cc
@@ -1,6 +1,6 @@
/*===========================================================================
- Copyright (C) 2006-2017 Yves Renard, Julien Pommier.
+ Copyright (C) 2006-2018 Yves Renard, Julien Pommier.
This file is a part of GetFEM++
@@ -36,12 +36,16 @@
#if GETFEM_HAVE_METIS_OLD_API
extern "C" void METIS_PartGraphKway(int *, int *, int *, int *, int *, int *,
int *, int *, int *, int *, int *);
-extern "C" void METIS_PartGraphRecursive(int *, int *, int *, int *, int *,
int *,
- int *, int *, int *, int *, int *);
-extern "C" void METIS_mCPartGraphKway(int *, int *, int *, int *, int *, int
*, int *,
- int *, int *, float *, int *, int *, int
*);
-extern "C" void METIS_mCPartGraphRecursive(int *, int *, int *, int *, int *,
int *, int *,
- int *, int *, int *, int *, int *);
+extern "C" void METIS_PartGraphRecursive(int *, int *, int *, int *, int *,
+ int *, int *, int *, int *, int *,
+ int *);
+extern "C" void METIS_mCPartGraphKway(int *, int *, int *, int *, int *, int *,
+ int *, int *, int *, float *, int *,
+ int *, int *);
+extern "C" void METIS_mCPartGraphRecursive(int *, int *, int *, int *, int *,
+ int *, int *, int *, int *, int *,
+ int *, int *);
+
#elif GETFEM_HAVE_METIS
# include <metis.h>
#endif
@@ -103,9 +107,8 @@ void asm_lsneuman_matrix
/**
- generic normal grad level set matrix (on the whole mesh level set or on
the specified
- convex set level set or boundary level set)
-
+ generic normal grad level set matrix (on the whole mesh level set or
+ on the specified element set level set or boundary level set)
*/
@@ -132,7 +135,7 @@ template<typename MAT> void asm_nlsgrad_matrix
/**************************************************************/
-/* assembling patch vector */
+/* assembling patch vector */
/**************************************************************/
template<class VEC>
@@ -150,7 +153,7 @@ void asm_patch_vector
}
/**************************************************************/
-/* assembling patch matrix */
+/* assembling patch matrix */
/**************************************************************/
template<class MAT>
@@ -430,6 +433,9 @@ static void do_high_level_generic_assembly(mexargs_in& in,
mexargs_out& out) {
bool with_secondary = false;
std::string secondary_domain;
+ bool with_select_output = false;
+ std::string select_var1, select_var2;
+
getfem::ga_workspace workspace2(md);
getfem::ga_workspace &workspace = with_model ? workspace2 : workspace1;
@@ -439,7 +445,14 @@ static void do_high_level_generic_assembly(mexargs_in& in,
mexargs_out& out) {
while (in.remaining()) {
std::string varname = in.pop().to_string();
- if (varname.compare("Secondary_domain") == 0 ||
+ if (varname.compare("select_output") == 0 ||
+ varname.compare("select output") == 0) {
+ GMM_ASSERT1(order > 0, "select_output option is for order 1 or 2"
+ "assemblies only");
+ with_select_output = true;
+ select_var1 = in.pop().to_string();
+ if (order == 2) select_var2 = in.pop().to_string();
+ } else if (varname.compare("Secondary_domain") == 0 ||
varname.compare("Secondary_Domain") == 0) {
GMM_ASSERT1(!with_secondary,
"Only one secondary domain can be specified");
@@ -502,7 +515,14 @@ static void do_high_level_generic_assembly(mexargs_in& in,
mexargs_out& out) {
getfem::model_real_plain_vector residual(nbdof);
workspace.set_assembled_vector(residual);
workspace.assembly(1);
- out.pop().from_dlvector(residual);
+ if (with_select_output) {
+ gmm::sub_interval I = workspace.interval_of_variable(select_var1);
+ getfem::model_real_plain_vector rresidual(I.size());
+ gmm::copy(gmm::sub_vector(residual, I), rresidual);
+ out.pop().from_dlvector(rresidual);
+ } else {
+ out.pop().from_dlvector(residual);
+ }
}
break;
@@ -511,9 +531,18 @@ static void do_high_level_generic_assembly(mexargs_in& in,
mexargs_out& out) {
getfem::model_real_sparse_matrix K(nbdof, nbdof);
workspace.set_assembled_matrix(K);
workspace.assembly(2);
- gf_real_sparse_by_col KK(nbdof, nbdof);
- gmm::copy(K, KK);
- out.pop().from_sparse(KK);
+
+ if (with_select_output) {
+ gmm::sub_interval I1 = workspace.interval_of_variable(select_var1);
+ gmm::sub_interval I2 = workspace.interval_of_variable(select_var2);
+ gf_real_sparse_by_col KK(I1.size(), I2.size());
+ gmm::copy(gmm::sub_matrix(K, I1, I2), KK);
+ out.pop().from_sparse(KK);
+ } else {
+ gf_real_sparse_by_col KK(nbdof, nbdof);
+ gmm::copy(K, KK);
+ out.pop().from_sparse(KK);
+ }
}
break;
@@ -738,7 +767,7 @@ void gf_asm(getfemint::mexargs_in& m_in,
getfemint::mexargs_out& m_out) {
if (subc_tab.size() == 0) {
- /address@hidden @CELL{...} = ('generic', @tmim mim, @int order, @str
expression, @int region, address@hidden model, ['Secondary_domain', 'name',]]
address@hidden varname, @int is_variable[, address@hidden mf, @tmimd mimd}],
value], ...)
+ /address@hidden @CELL{...} = ('generic', @tmim mim, @int order, @str
expression, @int region, address@hidden model, ['Secondary_domain', 'name',]]
address@hidden varname, @int is_variable[, address@hidden mf, @tmimd mimd}],
value], ['select_output', 'varname1'[, 'varname2]], ...)
High-level generic assembly procedure for volumic or boundary assembly.
Performs the generic assembly of `expression` with the integration
@@ -769,6 +798,12 @@ void gf_asm(getfemint::mexargs_in& m_in,
getfemint::mexargs_out& m_out) {
variables only (see GetFEM++ user documentation). Test functions are
only available for variables, not for constants.
+ `select_output` is an optional parameter which allows to reduce the
+ output vecotr (for `order` equal to 1) or the matrix (for `order`
+ equal to 2) to the degrees of freedom of the specified variables.
+ One variable has to be specified for a vector ouptut and two for a
+ matrix output.
+
Note that if several variables are given, the assembly of the
tangent matrix/residual vector will be done considering the order
in the call of the function (the degrees of freedom of the first
diff --git a/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
b/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
index af7e7f1..0e0cd0e 100644
--- a/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
+++ b/interface/tests/matlab/demo_Nitsche_contact_by_hand.m
@@ -43,19 +43,19 @@ ref_sol = 0 % 0 : Reference solution (Von Mises)
N = 2 % 2 or 3 dimensions
-R=0.25; % Radiaus of Omega_1.
+R=0.25; % Radius of Omega_1.
dirichlet_val = 0; % Dirchelet condition.
f_coeff=0; % friction coefficient.
clambda = 1; % Lame coefficient lambda.
cmu = 1; % Lame coefficient mu.
vertical_force = -0.1; % Vertical volume density of force on Omega_1.
penalty_parameter = 1E-7; % penalisation parmeter on Omega_1.
-elements_degree = 2 % degre of elments (1 or 2).
+elements_degree = 2 % degree of elements (1 or 2).
if (ref_sol == 0)
Theta = [-1]; % theta
Gamma0 = [1/100]; % Nitsche's parmeter gamma0
- Nxy =[50]; % mesh size (=1/nxy) 2D ->250 and 3D -> 25
+ Nxy =[50]; % mesh size (=1/nxy) 2D ->250 and 3D -> 25
elseif (ref_sol == 1)
diff --git a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index f07042c..2d7d6dc 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -1202,7 +1202,25 @@ namespace getfem {
: t(t_), inin(inin_), pt_type(ind_), nb(nb_) {}
};
-
+ struct ga_instruction_copy_interpolated_small_vect : public ga_instruction {
+ base_tensor &t;
+ const base_small_vector &vec;
+ ga_instruction_set::interpolate_info &inin;
+
+ virtual int exec() {
+ GA_DEBUG_INFO("Instruction: copy small vector");
+ GMM_ASSERT1(inin.ctx.is_convex_num_valid(), "Invalid element, "
+ "probably transformation failed");
+ GMM_ASSERT1(t.size() == vec.size(), "Invalid vector size.");
+ gmm::copy(vec, t.as_vector());
+ return 0;
+ }
+ ga_instruction_copy_interpolated_small_vect(base_tensor &t_,
+ const base_small_vector &vec_,
+ ga_instruction_set::interpolate_info &inin_)
+ : t(t_), vec(vec_), inin(inin_) {}
+ };
+
struct ga_instruction_interpolate : public ga_instruction {
base_tensor &t;
const mesh **m;
@@ -3789,7 +3807,6 @@ namespace getfem {
base_node P_ref;
size_type cv;
short_type face_num;
- gmm::clear(inin.Normal);
inin.pt_type = trans->transform(workspace, m, ctx, Normal, &(inin.m), cv,
face_num, P_ref, inin.Normal,
inin.derivatives, compute_der);
@@ -3807,6 +3824,7 @@ namespace getfem {
inin.pt_y = inin.ctx.xreal();
} else {
inin.ctx.invalid_convex_num();
+ inin.Normal.resize(0);
inin.pt_y = P_ref;
inin.has_ctx = false;
}
@@ -5010,9 +5028,10 @@ namespace getfem {
if (pnode->tensor().size() != m.dim())
pnode->init_vector_tensor(m.dim());
if (pnode->node_type == GA_NODE_INTERPOLATE_X)
- pgai = std::make_shared<ga_instruction_copy_small_vect>
+ pgai = std::make_shared<ga_instruction_copy_interpolated_small_vect>
(pnode->tensor(),
- rmi.interpolate_infos[pnode->interpolate_name].pt_y);
+ rmi.interpolate_infos[pnode->interpolate_name].pt_y,
+ rmi.interpolate_infos[pnode->interpolate_name]);
else if (pnode->node_type == GA_NODE_INTERPOLATE_NORMAL)
pgai = std::make_shared<ga_instruction_copy_Normal>
(pnode->tensor(),