branch: devel-logari81
commit 26f216c888fbd0edf13953fe2c6ad8012327f02a
Author: Konstantinos Poulios <logar...@gmail.com>
Date: Mon Aug 7 01:17:01 2017 +0200
white space and minor coding style normalization
---
src/bgeot_convex_ref.cc | 21 +-
src/bgeot_convex_structure.cc | 51 +-
src/bgeot_geometric_trans.cc | 14 +-
src/bgeot_poly.cc | 56 +-
src/getfem/bgeot_convex_ref.h | 16 +-
src/getfem/bgeot_convex_structure.h | 2 +-
src/getfem/bgeot_geometric_trans.h | 12 +-
src/getfem/bgeot_kdtree.h | 29 +-
src/getfem_contact_and_friction_common.cc | 8 +-
src/getfem_export.cc | 20 +-
src/getfem_fem.cc | 80 +--
src/getfem_import.cc | 32 +-
src/getfem_integration.cc | 869 +++++++++++++++---------------
src/getfem_mesh.cc | 4 +-
src/getfem_mesh_fem.cc | 118 ++--
15 files changed, 678 insertions(+), 654 deletions(-)
diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index 53c779f..2b3abe4 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -444,15 +444,18 @@ namespace bgeot {
dal::pstatic_stored_object_key
pk = std::make_shared<product_ref_key_>(a, b);
dal::pstatic_stored_object o = dal::search_stored_object(pk);
- if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
- pconvex_ref p = std::make_shared<product_ref_>(a, b);
- dal::add_stored_object(pk, p, a, b,
- convex_product_structure(a->structure(),
- b->structure()),
- p->pspt(), dal::PERMANENT_STATIC_OBJECT);
- pconvex_ref p1 = basic_convex_ref(p);
- if (p != p1) add_dependency(p, p1);
- return p;
+ if (o)
+ return std::dynamic_pointer_cast<const convex_of_reference>(o);
+ else {
+ pconvex_ref p = std::make_shared<product_ref_>(a, b);
+ dal::add_stored_object(pk, p, a, b,
+ convex_product_structure(a->structure(),
+ b->structure()),
+ p->pspt(), dal::PERMANENT_STATIC_OBJECT);
+ pconvex_ref p1 = basic_convex_ref(p);
+ if (p != p1) add_dependency(p, p1);
+ return p;
+ }
}
pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k) {
diff --git a/src/bgeot_convex_structure.cc b/src/bgeot_convex_structure.cc
index bb936f3..fe14175 100644
--- a/src/bgeot_convex_structure.cc
+++ b/src/bgeot_convex_structure.cc
@@ -299,7 +299,7 @@ namespace bgeot {
: convex_product_structure(cv1->faces_structure()[i], cv2);
faces[i] = std::vector<short_type>(cv1->nb_points_of_face(i)
- * cv2->nb_points());
+ * cv2->nb_points());
for (short_type j = 0; j < cv1->nb_points_of_face(i); j++)
for (short_type l = 0; l < cv2->nb_points(); l++) {
@@ -361,18 +361,23 @@ namespace bgeot {
DAL_DOUBLE_KEY(parallelepiped_key_, dim_type, dim_type);
pconvex_structure parallelepiped_structure(dim_type nc, dim_type k) {
- if (nc <= 1) return simplex_structure(nc, k);
+ if (nc <= 1)
+ return simplex_structure(nc, k);
+
dal::pstatic_stored_object_key
pcsk = std::make_shared<parallelepiped_key_>(nc, k);
dal::pstatic_stored_object o = dal::search_stored_object(pcsk);
- if (o) return ((std::dynamic_pointer_cast<const parallelepiped_>(o))->p);
- auto p = std::make_shared<parallelepiped_>();
- p->p = convex_product_structure(parallelepiped_structure(dim_type(nc-1),k),
- simplex_structure(1,k));
- dal::add_stored_object(pcsk, dal::pstatic_stored_object(p), p->p,
- dal::PERMANENT_STATIC_OBJECT);
- return p->p;
+ if (o)
+ return ((std::dynamic_pointer_cast<const parallelepiped_>(o))->p);
+ else {
+ auto p = std::make_shared<parallelepiped_>();
+ p->p =
convex_product_structure(parallelepiped_structure(dim_type(nc-1),k),
+ simplex_structure(1,k));
+ dal::add_stored_object(pcsk, dal::pstatic_stored_object(p), p->p,
+ dal::PERMANENT_STATIC_OBJECT);
+ return p->p;
+ }
}
/* ******************************************************************** */
@@ -398,7 +403,7 @@ namespace bgeot {
p->Nc = nc;
p->nbpt = (nc == 2) ? 8 : 20;
p->nbf = (nc == 2) ? 4 : 6;
- p->basic_pcvs = parallelepiped_structure(nc);
+ p->basic_pcvs = parallelepiped_structure(nc); // k=1
p->faces_struct.resize(p->nbf);
p->faces = std::vector< std::vector<short_type> >(p->nbf);
p->dir_points_ = std::vector<short_type>(p->Nc + 1);
@@ -468,24 +473,24 @@ namespace bgeot {
pconvex_structure pyramidal_structure(dim_type k) {
GMM_ASSERT1(k == 1 || k == 2, "Sorry, pyramidal elements implemented "
- "only for degree one or two.");
+ "only for degree one or two.");
dal::pstatic_stored_object_key
pcsk = std::make_shared<pyramidal_structure_key_>(k);
dal::pstatic_stored_object o = dal::search_stored_object(pcsk);
- if (o) return std::dynamic_pointer_cast<const convex_structure>(o);
+ if (o)
+ return std::dynamic_pointer_cast<const convex_structure>(o);
auto p = std::make_shared<pyramidal_structure_>();
pconvex_structure pcvs(p);
-
+
p->Nc = 3;
p->dir_points_ = std::vector<short_type>(p->Nc + 1);
-
if (k == 1) {
p->nbpt = 5;
p->nbf = 5;
p->auto_basic = true;
- // 4
+ // 4
// /|||
// / || |
// 2-|--|-3
@@ -508,12 +513,12 @@ namespace bgeot {
p->faces_struct[0] = parallelepiped_structure(2);
for (int i = 1; i < p->nbf; i++)
- p->faces_struct[i] = simplex_structure(2);
+ p->faces_struct[i] = simplex_structure(2);
dal::add_stored_object(pcsk, pcvs, parallelepiped_structure(2),
- simplex_structure(2),
- dal::PERMANENT_STATIC_OBJECT);
-
+ simplex_structure(2),
+ dal::PERMANENT_STATIC_OBJECT);
+
} else {
p->nbpt = 14;
p->nbf = 5;
@@ -544,15 +549,15 @@ namespace bgeot {
p->faces_struct[0] = parallelepiped_structure(2, 2);
for (int i = 1; i < p->nbf; i++)
- p->faces_struct[i] = simplex_structure(2, 2);
+ p->faces_struct[i] = simplex_structure(2, 2);
dal::add_stored_object(pcsk, pcvs, parallelepiped_structure(2, 2),
- simplex_structure(2, 2),
- dal::PERMANENT_STATIC_OBJECT);
+ simplex_structure(2, 2),
+ dal::PERMANENT_STATIC_OBJECT);
}
return pcvs;
}
-
+
/* ******************************************************************** */
/* Generic dummy convex with n global nodes. */
/* ******************************************************************** */
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 1d611e5..b07fc04 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -472,8 +472,8 @@ namespace bgeot {
{ vertex = false; break; }
if (vertex) vertices_.push_back(ip);
}
- auto dimension = dim();
- if (dynamic_cast<const torus_geom_trans *>(this)) dimension -= 1;
+ dim_type dimension = dim();
+ if (dynamic_cast<const torus_geom_trans *>(this)) --dimension;
GMM_ASSERT1(vertices_.size() > dimension, "Internal error");
}
@@ -801,7 +801,7 @@ namespace bgeot {
"2*x*x*y + 2*x*y*y - 3*x*y;");
for (int i = 0; i < 8; ++i)
- trans[i] = bgeot::read_base_poly(2, s);
+ trans[i] = read_base_poly(2, s);
} else {
std::stringstream s
("1 + 2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2"
@@ -831,7 +831,7 @@ namespace bgeot {
"2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2 - 5*x*y*z;");
for (int i = 0; i < 20; ++i)
- trans[i] = bgeot::read_base_poly(3, s);
+ trans[i] = read_base_poly(3, s);
}
fill_standard_vertices();
}
@@ -907,14 +907,14 @@ namespace bgeot {
};
static pgeometric_trans
- pyramidal_gt(gt_param_list& params,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ pyramidal_gt(gt_param_list& params,
+ std::vector<dal::pstatic_stored_object> &deps) {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
<< params.size() << " should be 1.");
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
int k = int(::floor(params[0].num() + 0.01));
- dependencies.push_back(pyramidal_element_of_reference(dim_type(k)));
+ deps.push_back(pyramidal_element_of_reference(dim_type(k)));
return std::make_shared<pyramidal_trans_>(dim_type(k));
}
diff --git a/src/bgeot_poly.cc b/src/bgeot_poly.cc
index de5ba06..545125e 100644
--- a/src/bgeot_poly.cc
+++ b/src/bgeot_poly.cc
@@ -33,9 +33,9 @@ namespace bgeot {
static bool init = false;
if (!init) {
for (short_type i = 0; i < STORED; ++i) {
- alpha_M_(i, 0) = alpha_M_(0, i) = 1;
- for (short_type j = 1; j <= i; ++j)
- alpha_M_(i,j) = alpha_M_(j,i) = (alpha_M_(i, j-1) * (i+j)) / j;
+ alpha_M_(i, 0) = alpha_M_(0, i) = 1;
+ for (short_type j = 1; j <= i; ++j)
+ alpha_M_(i,j) = alpha_M_(j,i) = (alpha_M_(i, j-1) * (i+j)) / j;
}
init = true;
}
@@ -46,7 +46,7 @@ namespace bgeot {
size_type alpha(short_type n, short_type d) {
alpha_init_();
GMM_ASSERT1(n < STORED && d < STORED,
- "alpha called with n = " << n << " and d = " << d);
+ "alpha called with n = " << n << " and d = " << d);
return alpha_(n,d);
}
@@ -56,7 +56,7 @@ namespace bgeot {
size_type g_idx = global_index_; short_type deg = degree_;
reverse_iterator it = rbegin() + 1;
for (l = short_type(n-2); l != short_type(-1); --l, ++it)
- if (*it != 0) break;
+ if (*it != 0) break;
short_type a = (*this)[n-1]; (*this)[n-1] = 0;
(*this)[short_type(l+1)] = short_type(a + 1);
if (l != short_type(-1)) ((*this)[l])--;
@@ -73,12 +73,12 @@ namespace bgeot {
size_type g_idx = global_index_; short_type deg = degree_;
reverse_iterator it = rbegin();
for (l = short_type(n-1); l != short_type(-1); --l, ++it) {
- if (*it != 0) break;
+ if (*it != 0) break;
}
if (l != short_type(-1)) {
- short_type a = (*this)[l];
- (*this)[l] = 0; (*this)[n-1] = short_type(a - 1);
- if (l > 0) ((*this)[l-1])++;
+ short_type a = (*this)[l];
+ (*this)[l] = 0; (*this)[n-1] = short_type(a - 1);
+ if (l > 0) ((*this)[l-1])++;
else if (short_type(deg+1)) degree_ = short_type(deg-1);
}
if (g_idx+1) global_index_ = g_idx-1;
@@ -144,20 +144,20 @@ namespace bgeot {
else if (s == "u" && n > 5) result = base_poly(n, 1, 5);
else if (s == "t" && n > 6) result = base_poly(n, 1, 6);
else if (s == "sqrt") {
- base_poly p = read_expression(n, f);
- if (p.degree() > 0) parse_error(1);
- result.one(); result *= sqrt(p[0]);
+ base_poly p = read_expression(n, f);
+ if (p.degree() > 0) parse_error(1);
+ result.one(); result *= sqrt(p[0]);
}
else { parse_error(2); }
break;
case 5 :
switch (s[0]) {
case '(' :
- result = read_base_poly(n, f);
- j = get_next_token(s, f);
- if (j != 5 || s[0] != ')') parse_error(3);
- break;
- default : parse_error(4);
+ result = read_base_poly(n, f);
+ j = get_next_token(s, f);
+ if (j != 5 || s[0] != ')') parse_error(3);
+ break;
+ default : parse_error(4);
}
break;
default : parse_error(5);
@@ -178,8 +178,8 @@ namespace bgeot {
}
void do_bin_op(std::vector<base_poly> &value_list,
- std::vector<int> &op_list,
- std::vector<int> &prior_list) {
+ std::vector<int> &op_list,
+ std::vector<int> &prior_list) {
base_poly &p2 = *(value_list.rbegin());
if (op_list.back() != 6) {
assert(value_list.size()>1);
@@ -190,14 +190,14 @@ namespace bgeot {
case 3 : p1 += p2; break;
case 4 : p1 -= p2; break;
case 5 :
- {
- if (p2.degree() > 0) parse_error(7);
- int pow = int(to_scalar(p2[0]));
- if (p2[0] != opt_long_scalar_type(pow) || pow < 0) parse_error(8);
- base_poly p = p1; p1.one();
- for (int i = 0; i < pow; ++i) p1 *= p;
- }
- break;
+ {
+ if (p2.degree() > 0) parse_error(7);
+ int pow = int(to_scalar(p2[0]));
+ if (p2[0] != opt_long_scalar_type(pow) || pow < 0) parse_error(8);
+ base_poly p = p1; p1.one();
+ for (int i = 0; i < pow; ++i) p1 *= p;
+ }
+ break;
default: assert(0);
}
value_list.pop_back();
@@ -223,7 +223,7 @@ namespace bgeot {
operator_priority_(i, i ? s[0] : '0', prior, op);
while (op) {
while (!prior_list.empty() && prior_list.back() <= prior)
- do_bin_op(value_list, op_list, prior_list);
+ do_bin_op(value_list, op_list, prior_list);
value_list.push_back(read_expression(n, f));
op_list.push_back(op);
diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h
index 8834afe..66db71e 100644
--- a/src/getfem/bgeot_convex_ref.h
+++ b/src/getfem/bgeot_convex_ref.h
@@ -111,12 +111,12 @@ namespace bgeot {
*/
virtual scalar_type is_in_face(short_type, const base_node &) const =0;
/// return the normal vector for each face.
- const std::vector<base_small_vector> &normals(void) const
+ const std::vector<base_small_vector> &normals() const
{ return normals_; }
/// return the vertices of the reference convex.
- const stored_point_tab &points(void) const { return *ppoints; }
- const stored_point_tab &points(void) { return *ppoints; }
- pstored_point_tab pspt(void) const { return ppoints; }
+ const stored_point_tab &points() const { return *ppoints; }
+ const stored_point_tab &points() { return *ppoints; }
+ pstored_point_tab pspt() const { return ppoints; }
virtual ~convex_of_reference()
{ DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Convex of reference"); }
@@ -130,7 +130,7 @@ namespace bgeot {
/// return the associated order 1 reference convex.
inline pconvex_ref basic_convex_ref(pconvex_ref cvr)
- { if (cvr->auto_basic) return cvr; else return cvr->basic_convex_ref_; }
+ { return cvr->auto_basic ? cvr : cvr->basic_convex_ref_; }
//@name public functions for obtaining a convex of reference
@@ -147,13 +147,11 @@ namespace bgeot {
*/
pconvex_ref Q2_incomplete_reference(dim_type d);
/** tensorial product of two convex ref.
- in order to ensure unicity, it is required the a->dim() >= b->dim()
- */
+ in order to ensure unicity, it is required the a->dim() >= b->dim() */
/** Pyramidal element of reference of degree k (k = 1 or 2 only) */
pconvex_ref pyramidal_element_of_reference(dim_type k);
pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b);
- /** equilateral simplex (degree 1). used only for mesh quality estimations
- */
+ /** equilateral simplex (degree 1). used only for mesh quality estimations */
pconvex_ref equilateral_simplex_of_reference(dim_type nc);
/** generic convex with n global nodes */
diff --git a/src/getfem/bgeot_convex_structure.h
b/src/getfem/bgeot_convex_structure.h
index 0dcb509..9a8cbe3 100644
--- a/src/getfem/bgeot_convex_structure.h
+++ b/src/getfem/bgeot_convex_structure.h
@@ -118,7 +118,7 @@ namespace bgeot {
/** Give a pointer array on the structures of the faces.
* faces_structure()[i] is a pointer on the structure of the face i.
*/
- inline const convex_structure_faces_ct &faces_structure(void) const
+ inline const convex_structure_faces_ct &faces_structure() const
{ return faces_struct; }
/** Return "direct" points indexes for a given face.
* @param i the face number.
diff --git a/src/getfem/bgeot_geometric_trans.h
b/src/getfem/bgeot_geometric_trans.h
index 01d1deb..9851024 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -126,7 +126,7 @@ namespace bgeot {
pconvex_structure structure() const { return cvr->structure(); }
/// Basic structure of the reference element.
pconvex_structure basic_structure() const
- { return bgeot::basic_structure(cvr->structure()); }
+ { return bgeot::basic_structure(cvr->structure()); }
/// Gives the value of the functions vector at a certain point.
virtual void poly_vector_val(const base_node &pt, base_vector &val) const
= 0;
/// Gives the value of a subgroup of the functions vector at a certain
point.
@@ -163,7 +163,7 @@ namespace bgeot {
if the transformation is linear, x is not used at all */
size_type complexity() const { return complexity_; }
virtual ~geometric_trans()
- { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Geometric transformation"); }
+ { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Geometric transformation"); }
geometric_trans()
{ DAL_STORED_OBJECT_DEBUG_CREATED(this, "Geometric transformation"); }
};
@@ -251,14 +251,14 @@ namespace bgeot {
*
* pt is the position of the evaluation point on the reference element
*/
- base_small_vector APIDECL compute_normal(const
geotrans_interpolation_context& c,
- size_type face);
+ base_small_vector APIDECL
+ compute_normal(const geotrans_interpolation_context& c, size_type face);
/** return the local basis (i.e. the normal in the first column, and the
* tangent vectors in the other columns
*/
- base_matrix
- APIDECL compute_local_basis(const geotrans_interpolation_context& c,
+ base_matrix APIDECL
+ compute_local_basis(const geotrans_interpolation_context& c,
size_type face);
//@}
diff --git a/src/getfem/bgeot_kdtree.h b/src/getfem/bgeot_kdtree.h
index 9d45555..81c00e9 100644
--- a/src/getfem/bgeot_kdtree.h
+++ b/src/getfem/bgeot_kdtree.h
@@ -52,7 +52,7 @@ namespace bgeot {
virtual ~kdtree_elt_base() {}
};
- /// store a point and the associated index for the kdtree.
+ /// store a point and the associated index for the kdtree.
/* std::pair<size_type,base_node> is not ok since it does not
have a suitable overloaded swap function ...
*/
@@ -69,15 +69,15 @@ namespace bgeot {
typedef std::vector<index_node_pair> kdtree_tab_type;
/** Balanced tree over a set of points.
-
+
Once the tree have been built, it is possible to query very
quickly for the list of points lying in a given box. Note that
this is not a dynamic structure: once you start to call
kdtree::points_in_box, you should not use anymore kdtree::add_point.
-
+
Here is an example of use (which tries to find the mapping between
the dof of the mesh_fem and the node numbers of its mesh):
-
+
@code
void dof_to_nodes(const getfem::mesh_fem &mf) {
const getfem::getfem_mesh &m = mf.linked_mesh();
@@ -94,9 +94,9 @@ namespace bgeot {
tree.points_in_box(t, P, P2);
if (t.size()) {
assert(t.size() == 1);
- cout << " dof " << d << " maps to mesh node " << t[0].i << "\n";
+ cout << " dof " << d << " maps to mesh node " << t[0].i << "\n";
}
- }
+ }
}
@endcode
*/
@@ -110,30 +110,31 @@ namespace bgeot {
void clear() { clear_tree(); pts = kdtree_tab_type(); N = 0; }
void reserve(size_type n) { pts.reserve(n); }
/// insert a new point
- size_type add_point(const base_node& n) {
+ size_type add_point(const base_node& n) {
size_type i = pts.size(); add_point_with_id(n,i); return i;
}
/// insert a new point, with an associated number.
void add_point_with_id(const base_node& n, size_type i) {
- if (pts.size() == 0) N = n.size();
- else if (N != n.size())
- GMM_ASSERT2(false, "invalid dimension");
+ if (pts.size() == 0)
+ N = n.size();
+ else
+ GMM_ASSERT2(N == n.size(), "invalid dimension");
if (tree) clear_tree();
pts.push_back(index_node_pair(i, n));
}
size_type nb_points() const { return pts.size(); }
const kdtree_tab_type &points() const { return pts; }
- /* fills ipts with the indexes of points in the box
+ /* fills ipts with the indexes of points in the box
[min,max] */
void points_in_box(kdtree_tab_type &ipts,
- const base_node &min,
- const base_node &max);
+ const base_node &min,
+ const base_node &max);
/* assigns at ipt the index of the nearest neighbor at location
pos and returns the square of the distance to this point*/
scalar_type nearest_neighbor(index_node_pair &ipt,
const base_node &pos);
private:
- typedef std::vector<size_type>::const_iterator ITER;
+ typedef std::vector<size_type>::const_iterator ITER;
void clear_tree();
};
}
diff --git a/src/getfem_contact_and_friction_common.cc
b/src/getfem_contact_and_friction_common.cc
index d65366d..385ce8c 100644
--- a/src/getfem_contact_and_friction_common.cc
+++ b/src/getfem_contact_and_friction_common.cc
@@ -1104,7 +1104,7 @@ namespace getfem {
for (size_type subiter(0);;) {
pps(a, hessa);
- det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1, false));
+ det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1,
false));
if (det > 1E-15) break;
for (size_type i = 0; i < N-1; ++i)
a[i] += gmm::random() * 1E-7;
@@ -1164,7 +1164,7 @@ namespace getfem {
for (size_type subiter(0);;) {
pps(a, hessa);
- det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1, false));
+ det = gmm::abs(bgeot::lu_inverse(&(*(hessa.begin())),N-1,
false));
if (det > 1E-15) break;
for (size_type i = 0; i < N-1; ++i)
a[i] += gmm::random() * 1E-7;
@@ -2677,8 +2677,8 @@ namespace getfem {
if (args.size() != 6) return false;
size_type N = args[0]->size();
if (N == 0 || args[1]->size() != N || args[2]->size() != N
- || args[3]->size() != 1 || args[4]->size() > 3
- || args[4]->size() == 0 || args[5]->size() != 1 ) return false;
+ || args[3]->size() != 1 || args[4]->size() > 3
+ || args[4]->size() == 0 || args[5]->size() != 1 ) return false;
ga_init_vector(sizes, N);
return true;
}
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index 14da099..7edd6cf 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -168,11 +168,11 @@ namespace getfem
pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt,
1)
: classical_fem(pgt, 1);
- short_type degree = 1;
- if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
- pgt->structure() != pgt->basic_structure())
- degree = 2;
-
+ short_type degree = 1;
+ if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
+ pgt->structure() != pgt->basic_structure())
+ degree = 2;
+
pmf->set_finite_element(cv, discontinuous ?
classical_discontinuous_fem(pgt, degree) :
classical_fem(pgt, degree));
@@ -480,9 +480,9 @@ namespace getfem
for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
bgeot::pgeometric_trans pgt2 = mf.linked_mesh().trans_of_convex(cv);
GMM_ASSERT1(basic_structure(pgt->structure()) ==
- basic_structure(pgt2->structure()),
- "Cannot export this mesh to opendx, it contains "
- "different convex types. Slice it first.");
+ basic_structure(pgt2->structure()),
+ "Cannot export this mesh to opendx, it contains "
+ "different convex types. Slice it first.");
pfem pf = mf.fem_of_element(cv);
bool discontinuous = false;
for (unsigned i=0; i < pf->nb_dof(cv); ++i) {
@@ -875,7 +875,7 @@ namespace getfem
if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
}
pfem classical_pf1 = discontinuous ?
- classical_discontinuous_fem(pgt, 1) : classical_fem(pgt, 1);
+ classical_discontinuous_fem(pgt, 1) : classical_fem(pgt, 1);
pmf->set_finite_element(cv, classical_pf1);
}
psl = NULL;
@@ -890,7 +890,7 @@ namespace getfem
case 4:
if ( 2 == pmf->fem_of_element(cv)->dim() ) t = POS_QU;
else if (3 == pmf->fem_of_element(cv)->dim()) t = POS_SI;
- break;
+ break;
case 6: t = POS_PR; break;
case 8: t = POS_HE; break;
case 5: t = POS_PY; break;
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 96a1b8f..5d81180 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -36,6 +36,8 @@
namespace getfem {
+ using bgeot::read_base_poly;
+ using bgeot::base_rational_fraction;
const base_matrix& fem_interpolation_context::M() const {
if (gmm::mat_nrows(M_) == 0) {
@@ -1157,7 +1159,7 @@ namespace getfem {
"4*(x*y - x*x*y);"
"2*x*x*y + 2*x*y*y - 3*x*y;");
- for (int i = 0; i < 8; ++i) p->base()[i] = bgeot::read_base_poly(2, s);
+ for (int i = 0; i < 8; ++i) p->base()[i] = read_base_poly(2, s);
p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.0));
p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.0));
@@ -1195,7 +1197,7 @@ namespace getfem {
"4*( - x^2*y*z + x*y*z);"
"2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2 - 5*x*y*z;");
- for (int i = 0; i < 20; ++i) p->base()[i] = bgeot::read_base_poly(3, s);
+ for (int i = 0; i < 20; ++i) p->base()[i] = read_base_poly(3, s);
p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.0, 0.0));
p->add_node(lagrange_dof(3), base_small_vector(0.5, 0.0, 0.0));
@@ -1256,7 +1258,7 @@ namespace getfem {
// 0---1---2
static pfem build_pyramidal_pk_fem(short_type k, bool disc) {
- auto p = std::make_shared<fem<bgeot::base_rational_fraction>>();
+ auto p = std::make_shared<fem<base_rational_fraction>>();
p->mref_convex() = bgeot::pyramidal_element_of_reference(1);
p->dim() = 3;
p->is_standard() = p->is_equivalent() = true;
@@ -1268,17 +1270,17 @@ namespace getfem {
if (k == 0) {
p->base().resize(1);
- p->base()[0] = bgeot::read_base_poly(3, "1");
+ p->base()[0] = read_base_poly(3, "1");
p->add_node(lagrange_0_dof(3), base_small_vector(0.0, 0.0, 0.5));
} else if (k == 1) {
p->base().resize(5);
- bgeot::base_rational_fraction // Q = xy/(1-z)
- Q(bgeot::read_base_poly(3, "x*y"), bgeot::read_base_poly(3, "1-z"));
- p->base()[0] = (bgeot::read_base_poly(3, "1-x-y-z") + Q)*0.25;
- p->base()[1] = (bgeot::read_base_poly(3, "1+x-y-z") - Q)*0.25;
- p->base()[2] = (bgeot::read_base_poly(3, "1-x+y-z") - Q)*0.25;
- p->base()[3] = (bgeot::read_base_poly(3, "1+x+y-z") + Q)*0.25;
- p->base()[4] = bgeot::read_base_poly(3, "z");
+ base_rational_fraction // Q = xy/(1-z)
+ Q(read_base_poly(3, "x*y"), read_base_poly(3, "1-z"));
+ p->base()[0] = (read_base_poly(3, "1-x-y-z") + Q)*0.25;
+ p->base()[1] = (read_base_poly(3, "1+x-y-z") - Q)*0.25;
+ p->base()[2] = (read_base_poly(3, "1-x+y-z") - Q)*0.25;
+ p->base()[3] = (read_base_poly(3, "1+x+y-z") + Q)*0.25;
+ p->base()[4] = read_base_poly(3, "z");
p->add_node(lag_dof, base_small_vector(-1.0, -1.0, 0.0));
p->add_node(lag_dof, base_small_vector( 1.0, -1.0, 0.0));
@@ -1289,16 +1291,16 @@ namespace getfem {
} else if (k == 2) {
p->base().resize(14);
- base_poly xi0 = bgeot::read_base_poly(3, "(1-z-x)*0.5");
- base_poly xi1 = bgeot::read_base_poly(3, "(1-z-y)*0.5");
- base_poly xi2 = bgeot::read_base_poly(3, "(1-z+x)*0.5");
- base_poly xi3 = bgeot::read_base_poly(3, "(1-z+y)*0.5");
- base_poly x = bgeot::read_base_poly(3, "x");
- base_poly y = bgeot::read_base_poly(3, "y");
- base_poly z = bgeot::read_base_poly(3, "z");
- base_poly ones = bgeot::read_base_poly(3, "1");
- base_poly un_z = bgeot::read_base_poly(3, "1-z");
- bgeot::base_rational_fraction Q(bgeot::read_base_poly(3, "1"), un_z);
+ base_poly xi0 = read_base_poly(3, "(1-z-x)*0.5");
+ base_poly xi1 = read_base_poly(3, "(1-z-y)*0.5");
+ base_poly xi2 = read_base_poly(3, "(1-z+x)*0.5");
+ base_poly xi3 = read_base_poly(3, "(1-z+y)*0.5");
+ base_poly x = read_base_poly(3, "x");
+ base_poly y = read_base_poly(3, "y");
+ base_poly z = read_base_poly(3, "z");
+ base_poly ones = read_base_poly(3, "1");
+ base_poly un_z = read_base_poly(3, "1-z");
+ base_rational_fraction Q(read_base_poly(3, "1"), un_z);
p->base()[ 0] = Q*Q*xi0*xi1*(x*y-z*un_z);
p->base()[ 1] = -Q*Q*xi0*xi1*xi2*y*4.;
@@ -1313,7 +1315,7 @@ namespace getfem {
p->base()[10] = Q*z*xi1*xi2*4.;
p->base()[11] = Q*z*xi3*xi0*4.;
p->base()[12] = Q*z*xi2*xi3*4.;
- p->base()[13] = bgeot::read_base_poly(3, "z*(2*z-1)");
+ p->base()[13] = read_base_poly(3, "z*(2*z-1)");
p->add_node(lag_dof, base_small_vector(-1.0, -1.0, 0.0));
p->add_node(lag_dof, base_small_vector( 0.0, -1.0, 0.0));
@@ -1338,8 +1340,7 @@ namespace getfem {
static pfem pyramidal_pk_fem
- (fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
short_type k = 2;
if (params.size() > 0) {
@@ -1347,14 +1348,13 @@ namespace getfem {
k = dim_type(::floor(params[0].num() + 0.01));
}
pfem p = build_pyramidal_pk_fem(k, false);
- dependencies.push_back(p->ref_convex(0));
- dependencies.push_back(p->node_tab(0));
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
return p;
}
static pfem pyramidal_disc_pk_fem
- (fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
short_type k = 2;
if (params.size() > 0) {
@@ -1362,18 +1362,18 @@ namespace getfem {
k = dim_type(::floor(params[0].num() + 0.01));
}
pfem p = build_pyramidal_pk_fem(k, true);
- dependencies.push_back(p->ref_convex(0));
- dependencies.push_back(p->node_tab(0));
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
return p;
}
- /* ******************************************************************** */
- /* P1 element with a bubble base fonction on a face */
- /* ******************************************************************** */
+ /* ******************************************************************** */
+ /* P1 element with a bubble base fonction on a face */
+ /* ******************************************************************** */
- struct P1_wabbfoaf_ : public PK_fem_ {
- P1_wabbfoaf_(dim_type nc);
- };
+ struct P1_wabbfoaf_ : public PK_fem_ {
+ P1_wabbfoaf_(dim_type nc);
+ };
P1_wabbfoaf_::P1_wabbfoaf_(dim_type nc) : PK_fem_(nc, 1) {
is_lag = false; es_degree = 2;
@@ -3136,7 +3136,7 @@ namespace getfem {
base_node pt(3);
for (unsigned k = 0; k < 5; ++k) {
for (unsigned i = 0; i < 4; ++i) {
- base_[k*4+i] = bgeot::read_base_poly(3, s);
+ base_[k*4+i] = read_base_poly(3, s);
pt[0] = pt[1] = pt[2] = ((k == 4) ? 1.0/3.0 : 0.0);
if (k == 4 && i) pt[i-1] = 0.0;
if (k < 4 && k) pt[k-1] = 1.0;
@@ -3308,7 +3308,7 @@ namespace getfem {
base_node pt(2);
for (unsigned k = 0; k < 7; ++k) {
for (unsigned i = 0; i < 3; ++i) {
- base_[k*3+i] = bgeot::read_base_poly(2, s);
+ base_[k*3+i] = read_base_poly(2, s);
if (k == 6) {
pt[0] = pt[1] = 0.5; if (i) pt[i-1] = 0.0;
add_node(normal_derivative_dof(2), pt);
@@ -3416,7 +3416,7 @@ namespace getfem {
"((x+y)^2 - x - y)*sqrt(2)/2; x*(x-1); y*(y-1);");
for (unsigned k = 0; k < 6; ++k)
- base_[k] = bgeot::read_base_poly(2, s);
+ base_[k] = read_base_poly(2, s);
add_node(lagrange_dof(2), base_node(0.0, 0.0));
add_node(lagrange_dof(2), base_node(1.0, 0.0));
@@ -3556,7 +3556,7 @@ namespace getfem {
std::stringstream name;
// Identifying if it is a torus structure
- if(bgeot::is_torus_geom_trans(pgt) && n == 3) n = 2;
+ if (bgeot::is_torus_geom_trans(pgt) && n == 3) n = 2;
/* Identifying P1-simplexes. */
if (nbp == n+1)
diff --git a/src/getfem_import.cc b/src/getfem_import.cc
index 6860d80..524f7ed 100644
--- a/src/getfem_import.cc
+++ b/src/getfem_import.cc
@@ -123,7 +123,7 @@ namespace getfem {
nodes.resize(6);
} break;
case 7: { /* PYRAMID */
- nodes.resize(5);
+ nodes.resize(5);
} break;
case 8: { /* 2ND ORDER LINE */
nodes.resize(3);
@@ -329,10 +329,10 @@ namespace getfem {
size_type j;
f >> j;
std::map<size_type, size_type>::iterator
- it = msh_node_2_getfem_node.find(j);
+ it = msh_node_2_getfem_node.find(j);
GMM_ASSERT1(it != msh_node_2_getfem_node.end(),
- "Invalid node ID " << j << " in gmsh element "
- << (ci.id + 1));
+ "Invalid node ID " << j << " in gmsh element "
+ << (ci.id + 1));
ci.nodes[i] = it->second;
}
if(ci.type != 15) ci.set_pgt();
@@ -715,11 +715,11 @@ namespace getfem {
for (size_type i=0; i < nnode; ++i) {
size_type j;
s >> j;
- std::map<size_type, size_type>::iterator
- it = msh_node_2_getfem_node.find(j);
+ std::map<size_type, size_type>::iterator
+ it = msh_node_2_getfem_node.find(j);
GMM_ASSERT1(it != msh_node_2_getfem_node.end(),
- "Invalid node ID " << j << " in GiD element " << cv_id);
- cv_nodes[i] = it->second;
+ "Invalid node ID " << j << " in GiD element " << cv_id);
+ cv_nodes[i] = it->second;
}
getfem_cv_nodes.resize(nnode);
for (size_type i=0; i < nnode; ++i) {
@@ -1105,7 +1105,7 @@ namespace getfem {
getfem_cv_nodes.begin()));
if (itype < elt_cnt.size())
elt_cnt[itype] += 1;
- }
+ }
}
}
}
@@ -1157,12 +1157,14 @@ namespace getfem {
/*NE: nombre d'elements (premier nombre du fichier .noboite)
NP: nombre de points (deuxieme nombre du fichier .noboite)
- ligne_debut_NP: la ligne commence la liste des coordonnees des points
dans le fichier .noboite*/
+ ligne_debut_NP: la ligne commence la liste des coordonnees des points
dans le
+ fichier .noboite*/
f >> NE>>NP;
ligne_debut_NP=NE*4/6+3;
- //passer 3 premiers lignes du fichier .noboite (la liste des elements
commence a la quatrieme ligne)
+ //passer 3 premiers lignes du fichier .noboite (la liste des elements
commence a la
+ //quatrieme ligne)
string contenu;
for (i=1; i<=ligne_debut_NP; i++){
getline(f, contenu);
@@ -1211,16 +1213,16 @@ namespace getfem {
if(fichier_GiD) // si l'ouverture a r�ussi
{
- // instructions
- fichier_GiD.close(); // on referme le fichier
+ // instructions
+ fichier_GiD.close(); // on referme le fichier
}
else // sinon
cerr << "Erreur � l'ouverture !" << endl;
if(f) // si l'ouverture a r�ussi
{
- // instructions
- //f.close(); // on referme le fichier
+ // instructions
+ //f.close(); // on referme le fichier
}
else // sinon
cerr << "Erreur � l'ouverture !" << endl;
diff --git a/src/getfem_integration.cc b/src/getfem_integration.cc
index 047fc3f..85f5233 100644
--- a/src/getfem_integration.cc
+++ b/src/getfem_integration.cc
@@ -33,10 +33,10 @@
namespace getfem {
/*
- * dummy integration method
+ * dummy integration method
*/
static pintegration_method im_none(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 0, "IM_NONE does not accept any parameter");
return std::make_shared<integration_method>();
}
@@ -49,7 +49,7 @@ namespace getfem {
hum->resize(i);
bgeot::power_index mi(P.dim()); mi[P.dim()-1] = P.degree();
for (size_type k = i; k > j; --k, --mi)
- (*hum)[k-1] = int_monomial(mi);
+ (*hum)[k-1] = int_monomial(mi);
}
base_poly::const_iterator it = P.begin(), ite = P.end();
std::vector<long_scalar_type>::const_iterator itb = int_monomials.begin();
@@ -68,11 +68,11 @@ namespace getfem {
hum->resize(i);
bgeot::power_index mi(P.dim()); mi[P.dim()-1] = P.degree();
for (size_type k = i; k > j; --k, --mi)
- (*hum)[k-1] = int_monomial_on_face(mi, f);
+ (*hum)[k-1] = int_monomial_on_face(mi, f);
}
base_poly::const_iterator it = P.begin(), ite = P.end();
std::vector<long_scalar_type>::const_iterator itb = hum->begin();
- for ( ; it != ite; ++it, ++itb)
+ for ( ; it != ite; ++it, ++itb)
res += long_scalar_type(*it) * long_scalar_type(*itb);
return res;
}
@@ -85,8 +85,8 @@ namespace getfem {
long_scalar_type int_monomial(const bgeot::power_index &power) const;
- long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
- short_type f) const;
+ long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
+ short_type f) const;
simplex_poly_integration_(bgeot::pconvex_structure c)
{ cvs = c; int_face_monomials.resize(c->nb_faces()); }
@@ -101,38 +101,39 @@ namespace getfem {
itme = power.end();
for ( ; itm != itme; ++itm)
for (int k = 1; k <= *itm; ++k, ++fa)
- res *= long_scalar_type(k) / long_scalar_type(fa);
-
+ res *= long_scalar_type(k) / long_scalar_type(fa);
+
for (int k = 0; k < cvs->dim(); k++) { res /= long_scalar_type(fa); fa++; }
return res;
}
-
+
long_scalar_type simplex_poly_integration_::int_monomial_on_face
(const bgeot::power_index &power, short_type f) const {
long_scalar_type res = LONG_SCAL(0);
-
+
if (f == 0 || power[f-1] == 0.0) {
res = (f == 0) ? sqrt(long_scalar_type(cvs->dim())) : LONG_SCAL(1);
short_type fa = 1;
bgeot::power_index::const_iterator itm = power.begin(),
- itme = power.end();
+ itme = power.end();
for ( ; itm != itme; ++itm)
- for (int k = 1; k <= *itm; ++k, ++fa)
- res *= long_scalar_type(k) / long_scalar_type(fa);
-
+ for (int k = 1; k <= *itm; ++k, ++fa)
+ res *= long_scalar_type(k) / long_scalar_type(fa);
+
for (int k = 1; k < cvs->dim(); k++) { res/=long_scalar_type(fa); fa++; }
}
return res;
}
- static pintegration_method exact_simplex(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ exact_simplex(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
- << params.size() << " should be 1.");
+ << params.size() << " should be 1.");
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 0 && n < 100 && double(n) == params[0].num(),
- "Bad parameters");
+ "Bad parameters");
dependencies.push_back(bgeot::simplex_structure(dim_type(n)));
ppoly_integration ppi = std::make_shared<simplex_poly_integration_>
(bgeot::simplex_structure(dim_type(n)));
@@ -148,8 +149,8 @@ namespace getfem {
long_scalar_type int_monomial(const bgeot::power_index &power) const;
- long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
- short_type f) const;
+ long_scalar_type int_monomial_on_face(const bgeot::power_index &power,
+ short_type f) const;
plyint_mul_structure_(ppoly_integration a, ppoly_integration b);
};
@@ -161,7 +162,7 @@ namespace getfem {
std::copy(power.begin() + cv1->dim(), power.end(), mi2.begin());
return cv1->int_monomial(mi1) * cv2->int_monomial(mi2);
}
-
+
long_scalar_type plyint_mul_structure_::int_monomial_on_face
(const bgeot::power_index &power, short_type f) const {
bgeot::power_index mi1(cv1->dim()), mi2(cv2->dim());
@@ -172,30 +173,31 @@ namespace getfem {
return cv1->int_monomial_on_face(mi1,f) * cv2->int_monomial(mi2);
else
return cv1->int_monomial(mi1)
- * cv2->int_monomial_on_face(mi2, short_type(f-nfx));
+ * cv2->int_monomial_on_face(mi2, short_type(f-nfx));
}
-
- plyint_mul_structure_::plyint_mul_structure_(ppoly_integration a,
- ppoly_integration b) {
+
+ plyint_mul_structure_::plyint_mul_structure_(ppoly_integration a,
+ ppoly_integration b) {
cv1 = a; cv2 = b;
cvs = bgeot::convex_product_structure(cv1->structure(),
- cv2->structure());
+ cv2->structure());
int_face_monomials.resize(cvs->nb_faces());
}
- static pintegration_method product_exact(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ product_exact(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
- GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
- "Bad type of parameters");
+ << params.size() << " should be 2.");
+ GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
+ "Bad type of parameters");
pintegration_method a = params[0].method();
pintegration_method b = params[1].method();
GMM_ASSERT1(a->type() == IM_EXACT && b->type() == IM_EXACT,
- "Bad parameters");
+ "Bad parameters");
dependencies.push_back(a); dependencies.push_back(b);
dependencies.push_back(bgeot::convex_product_structure(a->structure(),
- b->structure()));
+ b->structure()));
ppoly_integration ppi = std::make_shared<plyint_mul_structure_>
(a->exact_method(), b->exact_method());
return std::make_shared<integration_method>(ppi);
@@ -205,36 +207,38 @@ namespace getfem {
/* integration on parallelepiped. */
/* ******************************************************************** */
- static pintegration_method exact_parallelepiped(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ static pintegration_method
+ exact_parallelepiped(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
- << params.size() << " should be 1.");
+ << params.size() << " should be 1.");
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 0 && n < 100 && double(n) == params[0].num(),
- "Bad parameters");
+ "Bad parameters");
std::stringstream name;
if (n == 1)
name << "IM_EXACT_SIMPLEX(1)";
- else
+ else
name << "IM_PRODUCT(IM_EXACT_PARALLELEPIPED(" << n-1
- << "),IM_EXACT_SIMPLEX(1)))";
+ << "),IM_EXACT_SIMPLEX(1)))";
return int_method_descriptor(name.str());
}
- static pintegration_method exact_prism(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ static pintegration_method
+ exact_prism(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
- << params.size() << " should be 1.");
+ << params.size() << " should be 1.");
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
- "Bad parameters");
+ "Bad parameters");
std::stringstream name;
name << "IM_PRODUCT(IM_EXACT_SIMPLEX(" << n-1
- << "),IM_EXACT_SIMPLEX(1))";
+ << "),IM_EXACT_SIMPLEX(1))";
return int_method_descriptor(name.str());
}
@@ -243,116 +247,116 @@ namespace getfem {
/* ********************************************************************* */
void approx_integration::add_point(const base_node &pt,
- scalar_type w,short_type f) {
+ scalar_type w,short_type f) {
GMM_ASSERT1(!valid, "Impossible to modify a valid integration method.");
if (gmm::abs(w) > 1.0E-15) {
++f;
GMM_ASSERT1(f <= cvr->structure()->nb_faces(), "Wrong argument.");
size_type i = pt_to_store[f].search_node(pt);
if (i == size_type(-1)) {
- i = pt_to_store[f].add_node(pt);
- int_coeffs.resize(int_coeffs.size() + 1);
- for (size_type j = f; j <= cvr->structure()->nb_faces(); ++j)
- repartition[j] += 1;
- for (size_type j = int_coeffs.size(); j > repartition[f]; --j)
- int_coeffs[j-1] = int_coeffs[j-2];
- int_coeffs[repartition[f]-1] = 0.0;
+ i = pt_to_store[f].add_node(pt);
+ int_coeffs.resize(int_coeffs.size() + 1);
+ for (size_type j = f; j <= cvr->structure()->nb_faces(); ++j)
+ repartition[j] += 1;
+ for (size_type j = int_coeffs.size(); j > repartition[f]; --j)
+ int_coeffs[j-1] = int_coeffs[j-2];
+ int_coeffs[repartition[f]-1] = 0.0;
}
int_coeffs[((f == 0) ? 0 : repartition[f-1]) + i] += w;
}
}
void approx_integration::add_point_norepeat(const base_node &pt,
- scalar_type w,
- short_type f) {
+ scalar_type w,
+ short_type f) {
short_type f2 = f; ++f2;
if (pt_to_store[f2].search_node(pt) == size_type(-1)) add_point(pt,w,f);
}
void approx_integration::add_point_full_symmetric(base_node pt,
- scalar_type w) {
+ scalar_type w) {
dim_type n = cvr->structure()->dim();
dim_type k;
base_node pt2(n);
if (n+1 == cvr->structure()->nb_faces()) {
// simplices
- // for a point at (x,y) you have to consider the other points
+ // for a point at (x,y) you have to consider the other points
// at (y,x) (x,1-x-y) (1-x-y,x) (y,1-x-y) (1-x-y,y)
base_node pt3(n+1);
std::copy(pt.begin(), pt.end(), pt3.begin());
pt3[n] = 1; for (k = 0; k < n; ++k) pt3[n] -= pt[k];
std::vector<int> ind(n); std::fill(ind.begin(), ind.end(), 0);
- std::vector<bool> ind2(n+1);
+ std::vector<bool> ind2(n+1);
for (;;) {
-
- std::fill(ind2.begin(), ind2.end(), false);
- bool good = true;
- for (k = 0; k < n; ++k)
- if (ind2[ind[k]]) { good = false; break; } else ind2[ind[k]] = true;
- if (good) {
- for (k = 0; k < n; ++k) pt2[k] = pt3[ind[k]];
- add_point_norepeat(pt2, w);
- }
- ind[0]++; k = 0;
- while(ind[k] == n+1) { ind[k++] = 0; if (k == n) return; ind[k]++; }
+
+ std::fill(ind2.begin(), ind2.end(), false);
+ bool good = true;
+ for (k = 0; k < n; ++k)
+ if (ind2[ind[k]]) { good = false; break; } else ind2[ind[k]] = true;
+ if (good) {
+ for (k = 0; k < n; ++k) pt2[k] = pt3[ind[k]];
+ add_point_norepeat(pt2, w);
+ }
+ ind[0]++; k = 0;
+ while(ind[k] == n+1) { ind[k++] = 0; if (k == n) return; ind[k]++; }
}
}
else if (cvr->structure()->nb_points() == (size_type(1) << n)) {
// parallelepidedic
for (size_type i = 0; i < (size_type(1) << n); ++i) {
- for (k = 0; k < n; ++k)
- if (i & (size_type(1) << k)) pt2[k]=pt[k]; else pt2[k] = 1.0-pt[k];
- add_point_norepeat(pt2, w);
+ for (k = 0; k < n; ++k)
+ if (i & (size_type(1) << k)) pt2[k]=pt[k]; else pt2[k] = 1.0-pt[k];
+ add_point_norepeat(pt2, w);
}
}
else
GMM_ASSERT1(false, "Fully symmetric option is only valid for"
- "simplices and parallelepipedic elements");
+ "simplices and parallelepipedic elements");
}
void approx_integration::add_method_on_face(pintegration_method ppi,
- short_type f) {
+ short_type f) {
papprox_integration pai = ppi->approx_method();
GMM_ASSERT1(!valid, "Impossible to modify a valid integration method.");
GMM_ASSERT1(pai->structure() == structure()->faces_structure()[f],
- "structures missmatch");
+ "structures missmatch");
GMM_ASSERT1(ppi->type() == IM_APPROX, "Impossible with an exact method.");
dim_type N = pai->structure()->dim();
scalar_type det = 1.0;
base_node pt(N+1);
- std::vector<base_node> pts(N);
+ std::vector<base_node> pts(N);
for (size_type i = 0; i < N; ++i)
pts[i] = (cvr->dir_points_of_face(f))[i+1]
- - (cvr->dir_points_of_face(f))[0];
+ - (cvr->dir_points_of_face(f))[0];
if (N) {
base_matrix a(N+1, N), b(N, N), tmp(N, N);
for (dim_type i = 0; i < N+1; ++i)
- for (dim_type j = 0; j < N; ++j)
- a(i, j) = pts[j][i];
-
+ for (dim_type j = 0; j < N; ++j)
+ a(i, j) = pts[j][i];
+
gmm::mult(gmm::transposed(a), a, b);
det = ::sqrt(gmm::abs(gmm::lu_det(b)));
}
for (size_type i = 0; i < pai->nb_points_on_convex(); ++i) {
pt = (cvr->dir_points_of_face(f))[0];
for (dim_type j = 0; j < N; ++j)
- pt += pts[j] * ((*(pai->pintegration_points()))[i])[j];
+ pt += pts[j] * ((*(pai->pintegration_points()))[i])[j];
add_point(pt, pai->coeff(i) * det, f);
}
}
- void approx_integration::valid_method(void) {
+ void approx_integration::valid_method() {
std::vector<base_node> ptab(int_coeffs.size());
// std::vector<scalar_type> int_coeffs2(int_coeffs);
size_type i = 0;
for (short_type f = 0; f <= cvr->structure()->nb_faces(); ++f) {
// size_type j = i;
for (PT_TAB::const_iterator it = pt_to_store[f].begin();
- it != pt_to_store[f].end(); ++it /* , ++j */) {
- // int_coeffs[i] = int_coeffs2[j];
- ptab[i++] = *it;
+ it != pt_to_store[f].end(); ++it /* , ++j */) {
+ // int_coeffs[i] = int_coeffs2[j];
+ ptab[i++] = *it;
}
}
GMM_ASSERT1(i == int_coeffs.size(), "internal error.");
@@ -365,60 +369,61 @@ namespace getfem {
/* ********************************************************************* */
/* methods stored in getfem_im_list.h */
/* ********************************************************************* */
-
+
/// search a method in getfem_im_list.h
- static pintegration_method im_list_integration(std::string name,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ im_list_integration(std::string name,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
// cerr << "searching " << name << endl;
for (int i = 0; i < NB_IM; ++i)
if (!name.compare(im_desc_tab[i].method_name)) {
- bgeot::pgeometric_trans pgt
- = bgeot::geometric_trans_descriptor(im_desc_tab[i].geotrans_name);
- dim_type N = pgt->structure()->dim();
- base_node pt(N);
- auto pai = std::make_shared<approx_integration>(pgt->convex_ref());
- size_type fr = im_desc_tab[i].firstreal;
- for (size_type j = 0; j < im_desc_tab[i].nb_points; ++j) {
- for (dim_type k = 0; k < N; ++k)
- pt[k] = atof(im_desc_real[fr+j*(N+1)+k]);
- // pt[k] = LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+k]);
-
- switch (im_desc_node_type[im_desc_tab[i].firsttype + j]) {
- case 2: {
+ bgeot::pgeometric_trans pgt
+ = bgeot::geometric_trans_descriptor(im_desc_tab[i].geotrans_name);
+ dim_type N = pgt->structure()->dim();
+ base_node pt(N);
+ auto pai = std::make_shared<approx_integration>(pgt->convex_ref());
+ size_type fr = im_desc_tab[i].firstreal;
+ for (size_type j = 0; j < im_desc_tab[i].nb_points; ++j) {
+ for (dim_type k = 0; k < N; ++k)
+ pt[k] = atof(im_desc_real[fr+j*(N+1)+k]);
+ // pt[k] = LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+k]);
+
+ switch (im_desc_node_type[im_desc_tab[i].firsttype + j]) {
+ case 2: {
base_node pt2(pt.size());
for (bgeot::permutation p(pt.size()); !p.finished(); ++p) {
p.apply_to(pt,pt2);
- pai->add_point_full_symmetric(pt2,
- atof(im_desc_real[fr+j*(N+1)+N]));
- // pai->add_point_full_symmetric(pt2,
- // LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
+ pai->add_point_full_symmetric(pt2,
+ atof(im_desc_real[fr+j*(N+1)+N]));
+ // pai->add_point_full_symmetric(pt2,
+ // LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
}
- } break;
- case 1: {
- pai->add_point_full_symmetric(pt,atof(im_desc_real[fr+j*(N+1)+N]));
- // pai->add_point_full_symmetric(pt,
- // LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
- } break;
- case 0: {
- pai->add_point(pt, atof(im_desc_real[fr+j*(N+1)+N]));
- // pai->add_point(pt,LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
- } break;
- default: GMM_ASSERT1(false, "");
- }
- }
-
- for (short_type f = 0; N > 0 && f < pgt->structure()->nb_faces(); ++f)
- pai->add_method_on_face
- (int_method_descriptor
- (im_desc_face_meth[im_desc_tab[i].firstface + f]), f);
-
- pai->valid_method();
+ } break;
+ case 1: {
+ pai->add_point_full_symmetric(pt,atof(im_desc_real[fr+j*(N+1)+N]));
+ // pai->add_point_full_symmetric(pt,
+ // LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
+ } break;
+ case 0: {
+ pai->add_point(pt, atof(im_desc_real[fr+j*(N+1)+N]));
+ // pai->add_point(pt,LONG_SCALAR_ATOF(im_desc_real[fr+j*(N+1)+N]));
+ } break;
+ default: GMM_ASSERT1(false, "");
+ }
+ }
+
+ for (short_type f = 0; N > 0 && f < pgt->structure()->nb_faces(); ++f)
+ pai->add_method_on_face
+ (int_method_descriptor
+ (im_desc_face_meth[im_desc_tab[i].firstface + f]), f);
+
+ pai->valid_method();
// cerr << "finding " << name << endl;
- pintegration_method p(std::make_shared<integration_method>(pai));
- dependencies.push_back(p->approx_method()->ref_convex());
- dependencies.push_back(p->approx_method()->pintegration_points());
- return p;
+ pintegration_method p(std::make_shared<integration_method>(pai));
+ dependencies.push_back(p->approx_method()->ref_convex());
+ dependencies.push_back(p->approx_method()->pintegration_points());
+ return p;
}
return pintegration_method();
}
@@ -450,7 +455,7 @@ namespace getfem {
(base_poly(1,1,0) * polynomials[nb_lp-1]
* ((2.0 * long_scalar_type(nb_lp) - 1.0) / long_scalar_type(nb_lp)))
+ (polynomials[nb_lp-2]
- * ((1.0 - long_scalar_type(nb_lp)) / long_scalar_type(nb_lp)));
+ * ((1.0 - long_scalar_type(nb_lp)) / long_scalar_type(nb_lp)));
roots[nb_lp].resize(nb_lp);
roots[nb_lp][nb_lp/2] = 0.0;
long_scalar_type a = -1.0, b, c, d, e, cv, ev, ecart, ecart2;
@@ -474,7 +479,7 @@ namespace getfem {
roots[nb_lp][nb_lp-k-1] = -c;
a = b;
}
- }
+ }
}
};
@@ -484,30 +489,30 @@ namespace getfem {
gauss_approx_integration_::gauss_approx_integration_(short_type nbpt) {
GMM_ASSERT1(nbpt <= 32000, "too much points");
-
+
cvr = bgeot::simplex_of_reference(1);
std::vector<base_node> int_points(nbpt+2);
int_coeffs.resize(nbpt+2);
repartition.resize(3);
- repartition[0] = nbpt;
+ repartition[0] = nbpt;
repartition[1] = nbpt + 1;
- repartition[2] = nbpt + 2;
-
+ repartition[2] = nbpt + 2;
+
Legendre_polynomials Lp;
Lp.init(nbpt);
-
+
for (short_type i = 0; i < nbpt; ++i) {
int_points[i].resize(1);
long_scalar_type lr = Lp.roots[nbpt][i];
int_points[i][0] = 0.5 + 0.5 * bgeot::to_scalar(lr);
int_coeffs[i] = bgeot::to_scalar((1.0 - gmm::sqr(lr))
- / gmm::sqr( long_scalar_type(nbpt)
- * (Lp.polynomials[nbpt-1].eval(&lr))));
+ / gmm::sqr( long_scalar_type(nbpt)
+ * (Lp.polynomials[nbpt-1].eval(&lr))));
}
-
+
int_points[nbpt].resize(1);
int_points[nbpt][0] = 1.0; int_coeffs[nbpt] = 1.0;
-
+
int_points[nbpt+1].resize(1);
int_points[nbpt+1][0] = 0.0; int_coeffs[nbpt+1] = 1.0;
pint_points = bgeot::store_point_tab(int_points);
@@ -515,14 +520,15 @@ namespace getfem {
}
- static pintegration_method gauss1d(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ gauss1d(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
- << params.size() << " should be 1.");
+ << params.size() << " should be 1.");
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n >= 0 && n < 32000 && double(n) == params[0].num(),
- "Bad parameters");
+ "Bad parameters");
if (n & 1) {
std::stringstream name;
name << "IM_GAUSS1D(" << n-1 << ")";
@@ -530,7 +536,7 @@ namespace getfem {
}
else {
papprox_integration
- pai = std::make_shared<gauss_approx_integration_>(short_type(n/2 + 1));
+ pai = std::make_shared<gauss_approx_integration_>(short_type(n/2 + 1));
pintegration_method p = std::make_shared<integration_method>(pai);
dependencies.push_back(p->approx_method()->ref_convex());
dependencies.push_back(p->approx_method()->pintegration_points());
@@ -553,97 +559,98 @@ namespace getfem {
: approx_integration(bgeot::simplex_of_reference(nc)) {
size_type R = bgeot::alpha(nc,k);
- base_node c(nc);
+ base_node c(nc);
if (nc == 0) {
add_point(c, scalar_type(1));
}
else {
-
+
std::stringstream name;
name << "IM_EXACT_SIMPLEX(" << int(nc) << ")";
ppoly_integration ppi =
int_method_descriptor(name.str())->exact_method();
-
+
size_type sum = 0, l;
c.fill(scalar_type(0.0));
if (k == 0) c.fill(1.0 / scalar_type(nc+1));
-
+
gmm::dense_matrix<long_scalar_type> M(R, R);
std::vector<long_scalar_type> F(R), U(R);
std::vector<bgeot::power_index> base(R);
std::vector<base_node> nodes(R);
-
+
bgeot::power_index pi(nc);
-
+
for (size_type r = 0; r < R; ++r, ++pi) {
- base[r] = pi; nodes[r] = c;
- if (k != 0 && nc > 0) {
- l = 0; c[l] += 1.0 / scalar_type(k); sum++;
- while (sum > k) {
- sum -= int(floor(0.5+(c[l] * k)));
- c[l] = 0.0; l++; if (l == nc) break;
- c[l] += 1.0 / scalar_type(k); sum++;
- }
- }
+ base[r] = pi; nodes[r] = c;
+ if (k != 0 && nc > 0) {
+ l = 0; c[l] += 1.0 / scalar_type(k); sum++;
+ while (sum > k) {
+ sum -= int(floor(0.5+(c[l] * k)));
+ c[l] = 0.0; l++; if (l == nc) break;
+ c[l] += 1.0 / scalar_type(k); sum++;
+ }
+ }
}
// if (nc == 1) {
-// M = bgeot::vsmatrix<long_scalar_type>((R+1)/2, (R+1)/2);
-// U = F = bgeot::vsvector<long_scalar_type>((R+1)/2);
-// gmm::clear(M);
+// M = bgeot::vsmatrix<long_scalar_type>((R+1)/2, (R+1)/2);
+// U = F = bgeot::vsvector<long_scalar_type>((R+1)/2);
+// gmm::clear(M);
// }
for (size_type r = 0; r < R; ++r) {
-// if (nc == 1) {
-// if (r < (R+1)/2) {
-// F[r] = ppi->int_monomial(base[R-1-r]);
-// cout << "F[" << r << "] = " << F[r] << endl;
-// }
-// }
-// else {
- F[r] = ppi->int_monomial(base[r]);
- //cout << "F[" << r << "] = " << F[r] << endl;
-// }
- for (size_type q = 0; q < R; ++q) {
-// if (nc == 1) {
-// if (r < (R+1)/2) {
-// if (q < (R+1)/2)
-// M(r, q) += bgeot::eval_monomial(base[R-1-r], nodes[q].begin());
-// else
-// M(r, R-1-q) += bgeot::eval_monomial(base[R-1-r],
-// nodes[q].begin());
-// }
-// }
-// else
- M(r, q) = bgeot::eval_monomial(base[r], nodes[q].begin());
- }
+// if (nc == 1) {
+// if (r < (R+1)/2) {
+// F[r] = ppi->int_monomial(base[R-1-r]);
+// cout << "F[" << r << "] = " << F[r] << endl;
+// }
+// }
+// else {
+ F[r] = ppi->int_monomial(base[r]);
+ //cout << "F[" << r << "] = " << F[r] << endl;
+// }
+ for (size_type q = 0; q < R; ++q) {
+// if (nc == 1) {
+// if (r < (R+1)/2) {
+// if (q < (R+1)/2)
+// M(r, q) += bgeot::eval_monomial(base[R-1-r],
nodes[q].begin());
+// else
+// M(r, R-1-q) += bgeot::eval_monomial(base[R-1-r],
+// nodes[q].begin());
+// }
+// }
+// else
+ M(r, q) = bgeot::eval_monomial(base[r], nodes[q].begin());
+ }
}
-
+
gmm::lu_solve(M, U, F);
for (size_type r = 0; r < R; ++r)
- add_point(nodes[r], bgeot::to_scalar(U[r]));
-
+ add_point(nodes[r], bgeot::to_scalar(U[r]));
+
std::stringstream name2;
name2 << "IM_NC(" << int(nc-1) << "," << int(k) << ")";
for (short_type f = 0; f < structure()->nb_faces(); ++f)
- add_method_on_face(int_method_descriptor(name2.str()), f);
+ add_method_on_face(int_method_descriptor(name2.str()), f);
}
valid_method();
}
- static pintegration_method Newton_Cotes(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ Newton_Cotes(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
- "Bad type of parameters");
+ "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n >= 0 && n < 100 && k >= 0 && k <= 150 &&
- double(n) == params[0].num() && double(k) == params[1].num(),
- "Bad parameters");
+ double(n) == params[0].num() && double(k) == params[1].num(),
+ "Bad parameters");
papprox_integration
pai = std::make_shared<Newton_Cotes_approx_integration_>(dim_type(n),
- short_type(k));
+ short_type(k));
pintegration_method p = std::make_shared<integration_method>(pai);
dependencies.push_back(p->approx_method()->ref_convex());
dependencies.push_back(p->approx_method()->pintegration_points());
@@ -661,7 +668,7 @@ namespace getfem {
a_int_pro_integration::a_int_pro_integration(papprox_integration a,
- papprox_integration b) {
+ papprox_integration b) {
cvr = bgeot::convex_ref_product(a->ref_convex(), b->ref_convex());
size_type n1 = a->nb_points_on_convex();
size_type n2 = b->nb_points_on_convex();
@@ -672,12 +679,12 @@ namespace getfem {
repartition[0] = n1 * n2;
for (size_type i1 = 0; i1 < n1; ++i1)
for (size_type i2 = 0; i2 < n2; ++i2) {
- int_coeffs[i1 + i2 * n1] = a->coeff(i1) * b->coeff(i2);
- int_points[i1 + i2 * n1].resize(dim());
- std::copy(a->point(i1).begin(), a->point(i1).end(),
- int_points[i1 + i2 * n1].begin());
- std::copy(b->point(i2).begin(), b->point(i2).end(),
- int_points[i1 + i2 * n1].begin() + a->dim());
+ int_coeffs[i1 + i2 * n1] = a->coeff(i1) * b->coeff(i2);
+ int_points[i1 + i2 * n1].resize(dim());
+ std::copy(a->point(i1).begin(), a->point(i1).end(),
+ int_points[i1 + i2 * n1].begin());
+ std::copy(b->point(i2).begin(), b->point(i2).end(),
+ int_points[i1 + i2 * n1].begin() + a->dim());
}
short_type f = 0;
for (short_type f1 = 0; f1 < a->structure()->nb_faces(); ++f1, ++f) {
@@ -688,16 +695,16 @@ namespace getfem {
int_points.resize(repartition[f+1]);
int_coeffs.resize(repartition[f+1]);
for (size_type i1 = 0; i1 < n1; ++i1)
- for (size_type i2 = 0; i2 < n2; ++i2) {
- int_coeffs[w + i1 + i2 * n1] = a->coeff_on_face(f1, i1)
- * b->coeff(i2);
- int_points[w + i1 + i2 * n1].resize(dim());
- std::copy(a->point_on_face(f1, i1).begin(),
- a->point_on_face(f1, i1).end(),
- int_points[w + i1 + i2 * n1].begin());
- std::copy(b->point(i2).begin(), b->point(i2).end(),
- int_points[w + i1 + i2 * n1].begin() + a->dim());
- }
+ for (size_type i2 = 0; i2 < n2; ++i2) {
+ int_coeffs[w + i1 + i2 * n1] = a->coeff_on_face(f1, i1)
+ * b->coeff(i2);
+ int_points[w + i1 + i2 * n1].resize(dim());
+ std::copy(a->point_on_face(f1, i1).begin(),
+ a->point_on_face(f1, i1).end(),
+ int_points[w + i1 + i2 * n1].begin());
+ std::copy(b->point(i2).begin(), b->point(i2).end(),
+ int_points[w + i1 + i2 * n1].begin() + a->dim());
+ }
}
for (short_type f2 = 0; f2 < b->structure()->nb_faces(); ++f2, ++f) {
n1 = a->nb_points_on_convex();
@@ -707,46 +714,48 @@ namespace getfem {
int_points.resize(repartition[f+1]);
int_coeffs.resize(repartition[f+1]);
for (size_type i1 = 0; i1 < n1; ++i1)
- for (size_type i2 = 0; i2 < n2; ++i2) {
- int_coeffs[w + i1 + i2 * n1] = a->coeff(i1)
- * b->coeff_on_face(f2, i2);
- int_points[w + i1 + i2 * n1].resize(dim());
- std::copy(a->point(i1).begin(), a->point(i1).end(),
- int_points[w + i1 + i2 * n1].begin());
- std::copy(b->point_on_face(f2, i2).begin(),
- b->point_on_face(f2, i2).end(),
- int_points[w + i1 + i2 * n1].begin() + a->dim());
- }
+ for (size_type i2 = 0; i2 < n2; ++i2) {
+ int_coeffs[w + i1 + i2 * n1] = a->coeff(i1)
+ * b->coeff_on_face(f2, i2);
+ int_points[w + i1 + i2 * n1].resize(dim());
+ std::copy(a->point(i1).begin(), a->point(i1).end(),
+ int_points[w + i1 + i2 * n1].begin());
+ std::copy(b->point_on_face(f2, i2).begin(),
+ b->point_on_face(f2, i2).end(),
+ int_points[w + i1 + i2 * n1].begin() + a->dim());
+ }
}
pint_points = bgeot::store_point_tab(int_points);
valid = true;
}
- static pintegration_method product_approx(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ product_approx(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
- "Bad type of parameters");
+ "Bad type of parameters");
pintegration_method a = params[0].method();
pintegration_method b = params[1].method();
GMM_ASSERT1(a->type() == IM_APPROX && b->type() == IM_APPROX,
- "Bad parameters");
+ "Bad parameters");
papprox_integration
pai = std::make_shared<a_int_pro_integration>(a->approx_method(),
- b->approx_method());
+ b->approx_method());
pintegration_method p = std::make_shared<integration_method>(pai);
dependencies.push_back(p->approx_method()->ref_convex());
dependencies.push_back(p->approx_method()->pintegration_points());
return p;
}
- static pintegration_method product_which(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ product_which(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
- "Bad type of parameters");
+ "Bad type of parameters");
pintegration_method a = params[0].method();
pintegration_method b = params[1].method();
if (a->type() == IM_EXACT || b->type() == IM_EXACT)
@@ -759,42 +768,44 @@ namespace getfem {
/* integration on parallelepiped with Newton Cotes formulae */
/* ********************************************************************* */
- static pintegration_method Newton_Cotes_para(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ static pintegration_method
+ Newton_Cotes_para(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
- "Bad type of parameters");
+ "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
- double(n) == params[0].num() && double(k) == params[1].num(),
- "Bad parameters");
+ double(n) == params[0].num() && double(k) == params[1].num(),
+ "Bad parameters");
std::stringstream name;
if (n == 1)
name << "IM_NC(1," << k << ")";
- else
+ else
name << "IM_PRODUCT(IM_NC_PARALLELEPIPED(" << n-1 << "," << k
- << "),IM_NC(1," << k << "))";
+ << "),IM_NC(1," << k << "))";
return int_method_descriptor(name.str());
}
- static pintegration_method Newton_Cotes_prism(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ static pintegration_method
+ Newton_Cotes_prism(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
- "Bad type of parameters");
+ "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
- double(n) == params[0].num() && double(k) == params[1].num(),
- "Bad parameters");
+ double(n) == params[0].num() && double(k) == params[1].num(),
+ "Bad parameters");
std::stringstream name;
name << "IM_PRODUCT(IM_NC(" << n-1 << "," << k << "),IM_NC(1,"
- << k << "))";
+ << k << "))";
return int_method_descriptor(name.str());
}
@@ -802,24 +813,25 @@ namespace getfem {
/* integration on parallelepiped with Gauss formulae */
/* ********************************************************************* */
- static pintegration_method Gauss_paramul(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ static pintegration_method
+ Gauss_paramul(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
- "Bad type of parameters");
+ "Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
- double(n) == params[0].num() && double(k) == params[1].num(),
- "Bad parameters");
+ double(n) == params[0].num() && double(k) == params[1].num(),
+ "Bad parameters");
std::stringstream name;
if (n == 1)
name << "IM_GAUSS1D(" << k << ")";
- else
+ else
name << "IM_PRODUCT(IM_GAUSS_PARALLELEPIPED(" << n-1 << "," << k
- << "),IM_GAUSS1D(" << k << "))";
+ << "),IM_GAUSS1D(" << k << "))";
return int_method_descriptor(name.str());
}
@@ -828,8 +840,8 @@ namespace getfem {
/* ******************************************************************** */
struct quasi_polar_integration : public approx_integration {
- quasi_polar_integration(papprox_integration base_im,
- size_type ip1, size_type ip2=size_type(-1)) :
+ quasi_polar_integration(papprox_integration base_im,
+ size_type ip1, size_type ip2=size_type(-1)) :
approx_integration
((base_im->structure() == bgeot::parallelepiped_structure(3)) ?
bgeot::pyramidal_element_of_reference(1)
@@ -839,11 +851,11 @@ namespace getfem {
enum { SQUARE, PRISM, TETRA_CYL, PRISM2, PYRAMID } what;
if (N == 2) what = SQUARE;
else if (base_im->structure() == bgeot::prism_structure(3))
- what = (ip2 == size_type(-1) || ip1 == ip2) ? PRISM2 : PRISM;
+ what = (ip2 == size_type(-1) || ip1 == ip2) ? PRISM2 : PRISM;
else if (base_im->structure() == bgeot::simplex_structure(3))
- what = TETRA_CYL;
+ what = TETRA_CYL;
else if (base_im->structure() == bgeot::parallelepiped_structure(3))
- what = PYRAMID;
+ what = PYRAMID;
else GMM_ASSERT1(false, "Incoherent integration method");
// The first geometric transformation collapse a face of
@@ -855,8 +867,8 @@ namespace getfem {
bgeot::pgeometric_trans pgt2 = bgeot::simplex_geotrans(N, 1);
std::vector<base_node> nodes2(N+1);
if (what == PYRAMID) {
- pgt2 = bgeot::pyramidal_geotrans(1);
- nodes2.resize(5);
+ pgt2 = bgeot::pyramidal_geotrans(1);
+ nodes2.resize(5);
}
std::vector<size_type> other_nodes; // for the construction of node2
bgeot::pgeometric_trans pgt3 = bgeot::simplex_geotrans(N, 1);
@@ -864,162 +876,163 @@ namespace getfem {
switch (what) {
case SQUARE :
- nodes1[3] = nodes1[1];
- nodes2[ip1] = nodes1[1]; ip2 = ip1;
- other_nodes.push_back(0);
- other_nodes.push_back(2);
- break;
+ nodes1[3] = nodes1[1];
+ nodes2[ip1] = nodes1[1]; ip2 = ip1;
+ other_nodes.push_back(0);
+ other_nodes.push_back(2);
+ break;
case PRISM :
- nodes1[4] = nodes1[0]; nodes1[5] = nodes1[1];
- nodes2[ip1] = nodes1[0];
- nodes2[ip2] = nodes1[1];
- other_nodes.push_back(2);
- other_nodes.push_back(6);
- break;
+ nodes1[4] = nodes1[0]; nodes1[5] = nodes1[1];
+ nodes2[ip1] = nodes1[0];
+ nodes2[ip2] = nodes1[1];
+ other_nodes.push_back(2);
+ other_nodes.push_back(6);
+ break;
case TETRA_CYL :
- nodes3[0] = nodes1[0]; nodes3[1] = nodes1[1];
- nodes3[2] = nodes1[2]; nodes3[3] = nodes1[4];
- // nodes1[4] = nodes1[0]; nodes1[7] = base_node(1.0, 1.0, 2.0);
- nodes2[ip1] = nodes1[1]; ip2 = ip1;
- other_nodes.push_back(0);
- other_nodes.push_back(2);
- other_nodes.push_back(4);
- break;
+ nodes3[0] = nodes1[0]; nodes3[1] = nodes1[1];
+ nodes3[2] = nodes1[2]; nodes3[3] = nodes1[4];
+ // nodes1[4] = nodes1[0]; nodes1[7] = base_node(1.0, 1.0, 2.0);
+ nodes2[ip1] = nodes1[1]; ip2 = ip1;
+ other_nodes.push_back(0);
+ other_nodes.push_back(2);
+ other_nodes.push_back(4);
+ break;
case PRISM2 :
- nodes2[ip1] = nodes1[4];
- other_nodes.push_back(0);
- other_nodes.push_back(1);
- other_nodes.push_back(2);
- break;
+ nodes2[ip1] = nodes1[4];
+ other_nodes.push_back(0);
+ other_nodes.push_back(1);
+ other_nodes.push_back(2);
+ break;
case PYRAMID :
- ip2 = ip1 = 0;
- nodes1[0] = base_node(-1.,-1., 0.);
- nodes1[1] = base_node( 1.,-1., 0.);
- nodes1[2] = base_node(-1., 1., 0.);
- nodes1[3] = base_node( 1., 1., 0.);
- nodes1[4] = base_node( 0., 0., 1.);
- nodes1[5] = nodes1[6] = nodes1[7] = nodes1[4];
- nodes2[ip1] = nodes1[0];
- other_nodes.push_back(4);
- other_nodes.push_back(3);
- other_nodes.push_back(2);
- other_nodes.push_back(1);
+ ip2 = ip1 = 0;
+ nodes1[0] = base_node(-1.,-1., 0.);
+ nodes1[1] = base_node( 1.,-1., 0.);
+ nodes1[2] = base_node(-1., 1., 0.);
+ nodes1[3] = base_node( 1., 1., 0.);
+ nodes1[4] = base_node( 0., 0., 1.);
+ nodes1[5] = nodes1[6] = nodes1[7] = nodes1[4];
+ nodes2[ip1] = nodes1[0];
+ other_nodes.push_back(4);
+ other_nodes.push_back(3);
+ other_nodes.push_back(2);
+ other_nodes.push_back(1);
}
for (size_type i = 0; i <= nodes2.size()-1; ++i)
- if (i != ip1 && i != ip2) {
- GMM_ASSERT3(!other_nodes.empty(), "");
- nodes2[i] = nodes1[other_nodes.back()];
- other_nodes.pop_back();
- }
+ if (i != ip1 && i != ip2) {
+ GMM_ASSERT3(!other_nodes.empty(), "");
+ nodes2[i] = nodes1[other_nodes.back()];
+ other_nodes.pop_back();
+ }
- base_matrix G1, G2, G3;
+ base_matrix G1, G2, G3;
base_matrix K(N, N), grad(N, N), K3(N, N), K4(N, N);
base_node normal1(N), normal2(N);
bgeot::geotrans_inv_convex gic(nodes2, pgt2);
scalar_type J1, J2, J3(1), J4(1);
-
+
bgeot::vectors_to_base_matrix(G1, nodes1);
bgeot::vectors_to_base_matrix(G2, nodes2);
for (size_type nc = 0; nc < 2; ++nc) {
-
- if (what == TETRA_CYL) {
- if (nc == 1) nodes3[0] = nodes1[3];
- bgeot::vectors_to_base_matrix(G3, nodes3);
- pgt3->poly_vector_grad(nodes1[0], grad);
- gmm::mult(gmm::transposed(grad), gmm::transposed(G3), K3);
- J3 = gmm::abs(gmm::lu_inverse(K3)); /* = 1 */
- }
-
- for (size_type i=0; i < base_im->nb_points(); ++i) {
-
- gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
-
- size_type fp = size_type(-1); // Search the face number in the
- if (i >= base_im->nb_points_on_convex()) { // original element
- size_type ii = i - base_im->nb_points_on_convex();
- for (short_type ff=0; ff < base_im->structure()->nb_faces(); ++ff) {
- if (ii < base_im->nb_points_on_face(ff)) { fp = ff; break; }
- else ii -= base_im->nb_points_on_face(ff);
- }
- normal1 = base_im->ref_convex()->normals()[fp];
- }
-
- base_node P = base_im->point(i);
- if (what == TETRA_CYL) {
- P = pgt3->transform(P, nodes3);
- scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
- K4(0, 1) = - y / one_minus_z;
- K4(1, 1) = 1.0 - x / one_minus_z;
- K4(2, 1) = - x * y / gmm::sqr(one_minus_z);
- J4 = gmm::abs(gmm::lu_det(K4));
- P[1] *= (1.0 - x / one_minus_z);
- }
- if (what == PRISM2) {
- scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
- K4(0,0) = one_minus_z; K4(2,0) = -x;
- K4(1,1) = one_minus_z; K4(2,1) = -y;
- J4 = gmm::abs(gmm::lu_det(K4));
- P[0] *= one_minus_z;
- P[1] *= one_minus_z;
- }
-
- base_node P1 = pgt1->transform(P, nodes1), P2(N);
- pgt1->poly_vector_grad(P, grad);
- gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
- J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
-
- if (fp != size_type(-1) && J1 > 1E-10 && J4 > 1E-10) {
- if (what == TETRA_CYL) {
- gmm::mult(K3, normal1, normal2);
- normal1 = normal2;
- }
- gmm::lu_inverse(K4);
- gmm::lu_inverse(K);
- gmm::mult(K4, normal1, normal2);
- gmm::mult(K, normal2, normal1);
- normal2 = normal1;
- J1 *= gmm::vect_norm2(normal2);
- normal2 /= gmm::vect_norm2(normal2);
- }
- gic.invert(P1, P2);
- GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
- "Point not in the convex ref : " << P2);
-
- pgt2->poly_vector_grad(P2, grad);
- gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
- J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
-
- if (i < base_im->nb_points_on_convex())
- add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
- else if (J1 > 1E-10) {
- short_type f = short_type(-1);
- for (short_type ff = 0; ff <= N; ++ff)
- if (gmm::abs(pgt2->convex_ref()->is_in_face(ff, P2)) < 1E-8) {
- GMM_ASSERT1(f == short_type(-1),
- "An integration point is common to two faces");
- f = ff;
- }
- if (f != short_type(-1)) {
- gmm::mult(K, normal2, normal1);
- add_point(P2,base_im->coeff(i)*J1*gmm::vect_norm2(normal1)/J2,f);
- }
- // else { cout << "Point " << P2 << " eliminated" << endl; }
- }
- }
- if (what != TETRA_CYL) break;
+
+ if (what == TETRA_CYL) {
+ if (nc == 1) nodes3[0] = nodes1[3];
+ bgeot::vectors_to_base_matrix(G3, nodes3);
+ pgt3->poly_vector_grad(nodes1[0], grad);
+ gmm::mult(gmm::transposed(grad), gmm::transposed(G3), K3);
+ J3 = gmm::abs(gmm::lu_inverse(K3)); /* = 1 */
+ }
+
+ for (size_type i=0; i < base_im->nb_points(); ++i) {
+
+ gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
+
+ size_type fp = size_type(-1); // Search the face number in the
+ if (i >= base_im->nb_points_on_convex()) { // original element
+ size_type ii = i - base_im->nb_points_on_convex();
+ for (short_type ff=0; ff < base_im->structure()->nb_faces(); ++ff)
{
+ if (ii < base_im->nb_points_on_face(ff)) { fp = ff; break; }
+ else ii -= base_im->nb_points_on_face(ff);
+ }
+ normal1 = base_im->ref_convex()->normals()[fp];
+ }
+
+ base_node P = base_im->point(i);
+ if (what == TETRA_CYL) {
+ P = pgt3->transform(P, nodes3);
+ scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
+ K4(0, 1) = - y / one_minus_z;
+ K4(1, 1) = 1.0 - x / one_minus_z;
+ K4(2, 1) = - x * y / gmm::sqr(one_minus_z);
+ J4 = gmm::abs(gmm::lu_det(K4));
+ P[1] *= (1.0 - x / one_minus_z);
+ }
+ if (what == PRISM2) {
+ scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
+ K4(0,0) = one_minus_z; K4(2,0) = -x;
+ K4(1,1) = one_minus_z; K4(2,1) = -y;
+ J4 = gmm::abs(gmm::lu_det(K4));
+ P[0] *= one_minus_z;
+ P[1] *= one_minus_z;
+ }
+
+ base_node P1 = pgt1->transform(P, nodes1), P2(N);
+ pgt1->poly_vector_grad(P, grad);
+ gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
+ J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
+
+ if (fp != size_type(-1) && J1 > 1E-10 && J4 > 1E-10) {
+ if (what == TETRA_CYL) {
+ gmm::mult(K3, normal1, normal2);
+ normal1 = normal2;
+ }
+ gmm::lu_inverse(K4);
+ gmm::lu_inverse(K);
+ gmm::mult(K4, normal1, normal2);
+ gmm::mult(K, normal2, normal1);
+ normal2 = normal1;
+ J1 *= gmm::vect_norm2(normal2);
+ normal2 /= gmm::vect_norm2(normal2);
+ }
+ gic.invert(P1, P2);
+ GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
+ "Point not in the convex ref : " << P2);
+
+ pgt2->poly_vector_grad(P2, grad);
+ gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
+ J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
+
+ if (i < base_im->nb_points_on_convex())
+ add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
+ else if (J1 > 1E-10) {
+ short_type f = short_type(-1);
+ for (short_type ff = 0; ff <= N; ++ff)
+ if (gmm::abs(pgt2->convex_ref()->is_in_face(ff, P2)) < 1E-8) {
+ GMM_ASSERT1(f == short_type(-1),
+ "An integration point is common to two faces");
+ f = ff;
+ }
+ if (f != short_type(-1)) {
+ gmm::mult(K, normal2, normal1);
+ add_point(P2,base_im->coeff(i)*J1*gmm::vect_norm2(normal1)/J2,f);
+ }
+ // else { cout << "Point " << P2 << " eliminated" << endl; }
+ }
+ }
+ if (what != TETRA_CYL) break;
}
valid_method();
}
};
- static pintegration_method quasi_polar(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ quasi_polar(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() >= 1 && params[0].type() == 1,
- "Bad parameters for quasi polar integration: the first "
- "parameter should be an integration method");
+ "Bad parameters for quasi polar integration: the first "
+ "parameter should be an integration method");
pintegration_method a = params[0].method();
GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
int ip1 = 0, ip2 = 0;
@@ -1027,31 +1040,32 @@ namespace getfem {
GMM_ASSERT1(params.size() == 1, "Bad number of parameters");
} else {
GMM_ASSERT1(params.size() == 2 || params.size() == 3,
- "Bad number of parameters : " << params.size()
- << " should be 2 or 3.");
+ "Bad number of parameters : " << params.size()
+ << " should be 2 or 3.");
GMM_ASSERT1(params[1].type() == 0
- && params.back().type() == 0, "Bad type of parameters");
+ && params.back().type() == 0, "Bad type of parameters");
ip1 = int(::floor(params[1].num() + 0.01));
ip2 = int(::floor(params.back().num() + 0.01));
}
int N = a->approx_method()->dim();
GMM_ASSERT1(N >= 2 && N <= 3 && ip1 >= 0 && ip2 >= 0 && ip1 <= N
- && ip2 <= N, "Bad parameters");
+ && ip2 <= N, "Bad parameters");
papprox_integration
pai = std::make_shared<quasi_polar_integration>(a->approx_method(),
- ip1, ip2);
+ ip1, ip2);
pintegration_method p = std::make_shared<integration_method>(pai);
dependencies.push_back(p->approx_method()->ref_convex());
dependencies.push_back(p->approx_method()->pintegration_points());
return p;
}
- static pintegration_method pyramid(im_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pintegration_method
+ pyramid(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
GMM_ASSERT1(params.size() == 1 && params[0].type() == 1,
- "Bad parameters for pyramid integration: the first "
- "parameter should be an integration method");
+ "Bad parameters for pyramid integration: the first "
+ "parameter should be an integration method");
pintegration_method a = params[0].method();
GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
int N = a->approx_method()->dim();
@@ -1071,8 +1085,9 @@ namespace getfem {
/* Naming system */
/* ******************************************************************** */
- pintegration_method structured_composite_int_method(im_param_list &,
- std::vector<dal::pstatic_stored_object> &);
+ pintegration_method
+ structured_composite_int_method(im_param_list &,
+ std::vector<dal::pstatic_stored_object> &);
pintegration_method HCT_composite_int_method(im_param_list ¶ms,
std::vector<dal::pstatic_stored_object> &dependencies);
@@ -1109,7 +1124,7 @@ namespace getfem {
};
pintegration_method int_method_descriptor(std::string name,
- bool throw_if_not_found) {
+ bool throw_if_not_found) {
size_type i = 0;
return dal::singleton<im_naming_system>::instance().method
(name, i, throw_if_not_found);
@@ -1122,13 +1137,13 @@ namespace getfem {
// allows the add of an integration method.
void add_integration_name(std::string name,
- dal::naming_system<integration_method>::pfunction f) {
+ dal::naming_system<integration_method>::pfunction f) {
dal::singleton<im_naming_system>::instance().add_suffix(name, f);
}
/* Fonctions pour la ref. directe. */
-
+
pintegration_method exact_simplex_im(size_type n) {
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pintegration_method, pim, 0);
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(size_type, d, -2);
@@ -1186,20 +1201,20 @@ namespace getfem {
if (nbp == n+1)
if (cvs == bgeot::simplex_structure(dim_type(n)))
- { name << "IM_EXACT_SIMPLEX("; found = true; }
-
+ { name << "IM_EXACT_SIMPLEX("; found = true; }
+
/* Identifying Q1-parallelepiped. */
if (!found && nbp == (size_type(1) << n))
if (cvs == bgeot::parallelepiped_structure(dim_type(n)))
- { name << "IM_EXACT_PARALLELEPIPED("; found = true; }
+ { name << "IM_EXACT_PARALLELEPIPED("; found = true; }
/* Identifying Q1-prisms. */
-
+
if (!found && nbp == 2 * n)
if (cvs == bgeot::prism_structure(dim_type(n)))
- { name << "IM_EXACT_PRISM("; found = true; }
-
+ { name << "IM_EXACT_PRISM("; found = true; }
+
// To be completed
if (found) {
@@ -1208,7 +1223,7 @@ namespace getfem {
cvs_last = cvs;
return im_last;
}
-
+
GMM_ASSERT1(false, "This element is not taken into account. Contact us");
}
@@ -1247,7 +1262,7 @@ namespace getfem {
<< " of degree >= " << int(degree));
} else if (cvs->is_product(&a,&b) ||
(bgeot::basic_structure(cvs).get() &&
bgeot::basic_structure(cvs)->is_product(&a,&b))) {
- name << "IM_PRODUCT("
+ name << "IM_PRODUCT("
<< name_of_int_method(classical_approx_im_(a,degree)) << ","
<< name_of_int_method(classical_approx_im_(b,degree)) << ")";
} else GMM_ASSERT1(false, "unknown convex structure!");
@@ -1255,7 +1270,7 @@ namespace getfem {
}
pintegration_method classical_approx_im(bgeot::pgeometric_trans pgt,
- dim_type degree) {
+ dim_type degree) {
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(bgeot::pgeometric_trans, pgt_last,
0);
DEFINE_STATIC_THREAD_LOCAL(dim_type, degree_last);
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pintegration_method, im_last, 0);
@@ -1267,13 +1282,13 @@ namespace getfem {
return im_last;
}
- pintegration_method im_none(void) {
+ pintegration_method im_none() {
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pintegration_method,im_last,0);
if (!im_last) im_last = int_method_descriptor("IM_NONE");
return im_last;
}
- /* try to integrate all monomials up to order 'order' and return the
+ /* try to integrate all monomials up to order 'order' and return the
max. error */
scalar_type test_integration_error(papprox_integration pim, dim_type order) {
short_type dim = pim->dim();
@@ -1282,10 +1297,10 @@ namespace getfem {
for (bgeot::power_index idx(dim); idx.degree() <= order; ++idx) {
opt_long_scalar_type sum(0), realsum;
for (size_type i=0; i < pim->nb_points_on_convex(); ++i) {
- opt_long_scalar_type prod = pim->coeff(i);
- for (size_type d=0; d < dim; ++d)
- prod *= pow(opt_long_scalar_type(pim->point(i)[d]), idx[d]);
- sum += prod;
+ opt_long_scalar_type prod = pim->coeff(i);
+ for (size_type d=0; d < dim; ++d)
+ prod *= pow(opt_long_scalar_type(pim->point(i)[d]), idx[d]);
+ sum += prod;
}
realsum = exact->exact_method()->int_monomial(idx);
error = std::max(error, gmm::abs(realsum-sum));
@@ -1295,7 +1310,7 @@ namespace getfem {
papprox_integration get_approx_im_or_fail(pintegration_method pim) {
GMM_ASSERT1(pim->type() == IM_APPROX, "error estimate work only with "
- "approximate integration methods");
+ "approximate integration methods");
return pim->approx_method();
}
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index caeb31d..84a21d5 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -299,8 +299,8 @@ namespace getfem {
}
size_type mesh::add_pyramid(size_type a, size_type b,
- size_type c, size_type d, size_type e) {
- size_type ipt[5]; ipt[0] = a; ipt[1] = b; ipt[2] = c; ipt[3] = d; ipt[4] =
e;
+ size_type c, size_type d, size_type e) {
+ size_type ipt[5] = {a, b, c, d, e};
return add_convex(bgeot::pyramidal_geotrans(1), &(ipt[0]));
}
diff --git a/src/getfem_mesh_fem.cc b/src/getfem_mesh_fem.cc
index c3b4927..0816453 100644
--- a/src/getfem_mesh_fem.cc
+++ b/src/getfem_mesh_fem.cc
@@ -143,20 +143,20 @@ namespace getfem {
" and target_dim " << int(pf->target_dim()) << " of " <<
name_of_fem(pf));
-
+
if (cv == f_elems.size()) {
- f_elems.push_back(pf);
- fe_convex.add(cv);
- dof_enumeration_made = false;
+ f_elems.push_back(pf);
+ fe_convex.add(cv);
+ dof_enumeration_made = false;
touch(); v_num = act_counter();
} else {
- if (cv > f_elems.size()) f_elems.resize(cv+1);
- if (!fe_convex.is_in(cv) || f_elems[cv] != pf) {
- fe_convex.add(cv);
- f_elems[cv] = pf;
- dof_enumeration_made = false;
- touch(); v_num = act_counter();
- }
+ if (cv > f_elems.size()) f_elems.resize(cv+1);
+ if (!fe_convex.is_in(cv) || f_elems[cv] != pf) {
+ fe_convex.add(cv);
+ f_elems[cv] = pf;
+ dof_enumeration_made = false;
+ touch(); v_num = act_counter();
+ }
}
}
}
@@ -335,7 +335,7 @@ namespace getfem {
// Information for global dof
dal::bit_vector encountered_global_dof, processed_elt;
dal::dynamic_array<size_type> ind_global_dof;
-
+
// Auxilliary variables
std::vector<size_type> itab;
base_node P(linked_mesh().dim());
@@ -349,25 +349,25 @@ namespace getfem {
bgeot::pgeotrans_precomp pgp = 0;
for (dal::bv_visitor cv(linked_mesh().convex_index());
- !cv.finished(); ++cv) {
+ !cv.finished(); ++cv) {
if (fe_convex.is_in(cv)) {
- gmm::copy(linked_mesh().points_of_convex(cv)[0], bmin);
- gmm::copy(bmin, bmax);
- for (size_type i = 0; i < linked_mesh().nb_points_of_convex(cv); ++i) {
- const base_node &pt = linked_mesh().points_of_convex(cv)[i];
- for (size_type d = 1; d < bmin.size(); ++d) {
- bmin[d] = std::min(bmin[d], pt[d]);
- bmax[d] = std::max(bmax[d], pt[d]);
- }
- }
- elt_car_sizes[cv] = gmm::vect_dist2_sqr(bmin, bmax);
+ gmm::copy(linked_mesh().points_of_convex(cv)[0], bmin);
+ gmm::copy(bmin, bmax);
+ for (size_type i = 0; i < linked_mesh().nb_points_of_convex(cv); ++i) {
+ const base_node &pt = linked_mesh().points_of_convex(cv)[i];
+ for (size_type d = 1; d < bmin.size(); ++d) {
+ bmin[d] = std::min(bmin[d], pt[d]);
+ bmax[d] = std::max(bmax[d], pt[d]);
+ }
+ }
+ elt_car_sizes[cv] = gmm::vect_dist2_sqr(bmin, bmax);
}
}
dal::bit_vector cv_done;
-
+
for (dal::bv_visitor cv(linked_mesh().convex_index());
- !cv.finished(); ++cv) { // Loop on elements
+ !cv.finished(); ++cv) { // Loop on elements
if (!fe_convex.is_in(cv)) continue;
pfem pf = fem_of_element(cv);
if (pf != first_pf) is_uniform_ = false;
@@ -400,34 +400,34 @@ namespace getfem {
pgp->transform(linked_mesh().points_of_convex(cv), i, P);
size_type idof = nbdof;
- if (dof_nodes[cv].nb_points() > 0) {
- scalar_type dist = dof_nodes[cv].nearest_neighbor(ipt, P);
- if (gmm::abs(dist) <= 1e-6*elt_car_sizes[cv]) {
- fd.ind_node=ipt.i;
- auto it = dof_sorts[cv].find(fd);
- if (it != dof_sorts[cv].end()) idof = it->second;
- }
- }
-
+ if (dof_nodes[cv].nb_points() > 0) {
+ scalar_type dist = dof_nodes[cv].nearest_neighbor(ipt, P);
+ if (gmm::abs(dist) <= 1e-6*elt_car_sizes[cv]) {
+ fd.ind_node=ipt.i;
+ auto it = dof_sorts[cv].find(fd);
+ if (it != dof_sorts[cv].end()) idof = it->second;
+ }
+ }
+
if (idof == nbdof) {
- nbdof += Qdim / pf->target_dim();
-
- linked_mesh().neighbours_of_convex(cv, pf->faces_of_dof(cv, i), s);
- for (size_type ncv : s) { // For each unscanned neighbour
- if (!cv_done[ncv] && fe_convex.is_in(ncv)) { // add the dof
-
- fd.ind_node = size_type(-1);
- if (dof_nodes[ncv].nb_points() > 0) {
- scalar_type dist = dof_nodes[ncv].nearest_neighbor(ipt, P);
- if (gmm::abs(dist) <= 1e-6*elt_car_sizes[ncv])
- fd.ind_node=ipt.i;
- }
- if (fd.ind_node == size_type(-1))
- fd.ind_node = dof_nodes[ncv].add_point(P);
- dof_sorts[ncv][fd] = idof;
- }
- }
- }
+ nbdof += Qdim / pf->target_dim();
+
+ linked_mesh().neighbours_of_convex(cv, pf->faces_of_dof(cv, i), s);
+ for (size_type ncv : s) { // For each unscanned neighbour
+ if (!cv_done[ncv] && fe_convex.is_in(ncv)) { // add the dof
+
+ fd.ind_node = size_type(-1);
+ if (dof_nodes[ncv].nb_points() > 0) {
+ scalar_type dist = dof_nodes[ncv].nearest_neighbor(ipt, P);
+ if (gmm::abs(dist) <= 1e-6*elt_car_sizes[ncv])
+ fd.ind_node=ipt.i;
+ }
+ if (fd.ind_node == size_type(-1))
+ fd.ind_node = dof_nodes[ncv].add_point(P);
+ dof_sorts[ncv][fd] = idof;
+ }
+ }
+ }
itab[i] = idof;
}
}
@@ -571,8 +571,8 @@ namespace getfem {
} else if (bgeot::casecmp(tmp, "DOF_ENUMERATION") == 0) {
dal::bit_vector doflst;
dof_structure.clear(); dof_enumeration_made = false;
- is_uniform_ = true;
- size_type nbdof_unif = size_type(-1);
+ is_uniform_ = true;
+ size_type nbdof_unif = size_type(-1);
touch(); v_num = act_counter();
while (true) {
bgeot::get_token(ist, tmp);
@@ -585,11 +585,11 @@ namespace getfem {
std::vector<size_type> tab;
if (convex_index().is_in(ic) && tmp.size() &&
isdigit(tmp[0]) && tmp2 == ":") {
- size_type nbd = nb_basic_dof_of_element(ic);
- if (nbdof_unif == size_type(-1))
- nbdof_unif = nbd;
- else if (nbdof_unif != nbd)
- is_uniform_ = false;
+ size_type nbd = nb_basic_dof_of_element(ic);
+ if (nbdof_unif == size_type(-1))
+ nbdof_unif = nbd;
+ else if (nbdof_unif != nbd)
+ is_uniform_ = false;
tab.resize(nb_basic_dof_of_element(ic));
for (size_type i=0; i < fem_of_element(ic)->nb_dof(ic);
i++) {
