GetFEM  5.4.2
getfem_torus.cc
1 /*===========================================================================
2 
3  Copyright (C) 2014-2020 Liang Jin Lim
4 
5  This file is a part of GetFEM
6 
7  GetFEM is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include "getfem/bgeot_torus.h"
23 #include "getfem/getfem_torus.h"
24 
25 namespace getfem
26 {
27 
28  void torus_fem::init(){
29  cvr = poriginal_fem_->ref_convex(0);
30  dim_ = cvr->structure()->dim();
31  is_standard_fem = false;
32  is_equiv = real_element_defined = true;
33  is_pol = poriginal_fem_->is_polynomial();
34  is_polycomp = poriginal_fem_->is_polynomialcomp();
35  is_lag = poriginal_fem_->is_lagrange();
36  es_degree = poriginal_fem_->estimated_degree();
37  ntarget_dim = 3;
38 
39  std::stringstream nm;
40  nm << "FEM_TORUS(" << poriginal_fem_->debug_name() << ")";
41  debug_name_ = nm.str();
42 
43  init_cvs_node();
44  GMM_ASSERT1(poriginal_fem_->target_dim() == 1, "Vectorial fems not supported");
45  size_type nb_dof_origin = poriginal_fem_->nb_dof(0);
46  for (size_type k = 0; k < nb_dof_origin; ++k)
47  {
48  for(size_type j = 0; j < 2; ++j){
49  add_node(xfem_dof(poriginal_fem_->dof_types()[k], j),
50  poriginal_fem_->node_of_dof(0, k));
51  }
52  }
53  }
54 
55  size_type torus_fem::index_of_global_dof
56  (size_type /*cv*/, size_type i) const
57  { return i; }
58 
59  void torus_fem::base_value(const base_node &, base_tensor &) const
60  { GMM_ASSERT1(false, "No base values, real only element."); }
61 
62  void torus_fem::grad_base_value(const base_node &, base_tensor &) const
63  { GMM_ASSERT1(false, "No grad values, real only element."); }
64 
65  void torus_fem::hess_base_value(const base_node &, base_tensor &) const{
66  GMM_ASSERT1(false, "No hess values, real only element.");
67  }
68 
70  base_tensor &t, bool) const{
71 
72  GMM_ASSERT1(!(poriginal_fem_->is_on_real_element()), "Original FEM must not be real.");
73 
74  base_tensor u_orig;
75  poriginal_fem_->base_value(c.xref(), u_orig);
76  if (!(poriginal_fem_->is_equivalent())){
77  base_tensor u_temp = u_orig;
78  u_orig.mat_transp_reduction(u_temp, c.M(), 0);
79  }
80 
81  if(is_scalar_){
82  t = u_orig;
83  return;
84  }
85 
86  //expand original base of [nb_base, 1]
87  //to vectorial form [nb_base * dim_, dim_ + 1]
88  bgeot::multi_index tensor_size(u_orig.sizes());
89  tensor_size[0] *= dim_;
90  tensor_size[1] = ntarget_dim;
91  t.adjust_sizes(tensor_size);
92  for (size_type i = 0; i < u_orig.sizes()[0]; ++i) {
93  for (dim_type j = 0; j < dim_; ++j) {
94  t(i*dim_ + j, j) = u_orig(i, 0);
95  }
96  }
97  }
98 
100  (const getfem::fem_interpolation_context& c, base_tensor &t, bool) const
101  {
102  GMM_ASSERT1(!(poriginal_fem_->is_on_real_element()), "Original FEM must not be real.");
103 
104  bgeot::scalar_type radius = std::abs(c.xreal()[0]);
105  base_tensor u;
106  bgeot::pstored_point_tab ppt = c.pgp()->get_ppoint_tab();
107  getfem::pfem_precomp pfp = getfem::fem_precomp(poriginal_fem_, ppt, 0);
108  base_tensor u_origin = pfp->grad(c.ii());
109  //poriginal_fem_->grad_base_value(c.xref(), u_origin);
110  GMM_ASSERT1(!u_origin.empty(), "Original FEM is unable to provide grad base value!");
111 
112  base_tensor n_origin = pfp->val(c.ii());
113  //poriginal_fem_->base_value(c.xref(), n_origin);
114  GMM_ASSERT1(!n_origin.empty(), "Original FEM is unable to provide base value!");
115 
116  //expand original grad of [nb_base, 1, dim_]
117  //to vectorial form [nb_base * dim_, dim_ + 1, dim_ + 1]
118  const bgeot::multi_index &origin_size = u_origin.sizes();
119  bgeot::multi_index tensor_size(origin_size);
120  dim_type dim_size = is_scalar_ ? 1 : dim_;
121  tensor_size[0] *= dim_size;
122  tensor_size[1] = ntarget_dim;
123  tensor_size[2] = dim_ + 1;
124  u.adjust_sizes(tensor_size);
125  for (size_type i = 0; i < origin_size[0]; ++i) { //dof
126  for (dim_type j = 0; j < dim_size; ++j) {
127  for (dim_type k = 0; k < dim_; ++k) {
128  u(i*dim_size+j, j, k) = u_origin(i, 0, k);
129  }
130  }
131  }
132  t = u;
133  t.mat_transp_reduction(u, c.B(), 2);
134 
135  if(is_scalar_) return;
136 
137  for (size_type i = 0; i < origin_size[0]; ++i) {
138  t(i*dim_size, dim_, dim_) = n_origin[i] / radius;
139  }
140  }
141 
143  (const getfem::fem_interpolation_context & /*c*/, base_tensor & /*t*/, bool) const
144  {
145  GMM_ASSERT1(false, "Hessian not yet implemented in torus fem.");
146  }
147 
148  void torus_fem::set_to_scalar(bool is_scalar){
149  if(is_scalar_ == is_scalar) return;
150 
151  is_scalar_ = is_scalar;
152 
153  if(is_scalar_){
154  ntarget_dim = 1;
155  dof_types_.clear();
156  init_cvs_node();
157  size_type nb_dof_origin = poriginal_fem_->nb_dof(0);
158  for (size_type k = 0; k < nb_dof_origin; ++k){
159  add_node(poriginal_fem_->dof_types()[k], poriginal_fem_->node_of_dof(0, k));
160  }
161  }else{
162  ntarget_dim = 3;
163  dof_types_.clear();
164  init_cvs_node();
165  size_type nb_dof_origin = poriginal_fem_->nb_dof(0);
166  for (size_type k = 0; k < nb_dof_origin; ++k)
167  {
168  for(size_type j = 0; j < 2; ++j){
169  add_node(xfem_dof(poriginal_fem_->dof_types()[k], j),
170  poriginal_fem_->node_of_dof(0, k));
171  }
172  }
173  }
174  }
175 
176  pfem torus_fem::get_original_pfem() const{return poriginal_fem_;}
177 
178  DAL_SIMPLE_KEY(torus_fem_key, bgeot::size_type);
179 
180  getfem::pfem new_torus_fem(getfem::pfem pf) {
181  static bgeot::size_type key_count = 0;
182  ++key_count;
183  getfem::pfem pfem_torus = std::make_shared<torus_fem>(pf);
184  dal::pstatic_stored_object_key
185  pk = std::make_shared<torus_fem_key>(key_count);
186  dal::add_stored_object(pk, pfem_torus, pfem_torus->node_tab(0));
187  return pfem_torus;
188  }
189 
190  void del_torus_fem(getfem::pfem pf){
191  const torus_fem *ptorus_fem = dynamic_cast<const torus_fem*>(pf.get());
192  if (ptorus_fem != 0) dal::del_stored_object(pf);
193  }
194 
195  void torus_mesh_fem::adapt_to_torus_(){
196 
197  for (dal::bv_visitor cv(linked_mesh().convex_index()); !cv.finished(); ++cv){
198  pfem poriginal_fem = fem_of_element(cv);
199  if(poriginal_fem == 0) continue;
200 
201  del_torus_fem(poriginal_fem);
202 
203  pfem pf = new_torus_fem(poriginal_fem);
204  torus_fem *pf_torus = dynamic_cast<torus_fem*>(const_cast<virtual_fem*>(pf.get()));
205  pf_torus->set_to_scalar((Qdim != 3));
206  set_finite_element(cv, pf);
207  }
208  touch();
209  }
210 
212  {
213  const_cast<torus_mesh_fem*>(this)->adapt_to_torus_();
214 
215  for (dal::bv_visitor cv(linked_mesh().convex_index()); !cv.finished(); ++cv){
216  pfem pf = fem_of_element(cv);
217  if(pf == 0) continue;
218  torus_fem *pf_torus = dynamic_cast<torus_fem*>(const_cast<virtual_fem*>(pf.get()));
219  if(pf_torus == 0) continue;
220  pf_torus->set_to_scalar((Qdim != 3));
221  }
222 
224  }
225 
226  torus_mesh::torus_mesh(std::string name) : mesh(std::move(name)){}
227 
229  {
230  base_matrix G;
231  bgeot::vectors_to_base_matrix(G, points_of_convex(ic));
232  G.resize(2, G.ncols());
233  auto pgt_torus = std::dynamic_pointer_cast<const bgeot::torus_geom_trans>(trans_of_convex(ic));
234  GMM_ASSERT2(pgt_torus, "Internal error, convex is not a torus transformation.");
235  return getfem::convex_radius_estimate(pgt_torus->get_original_transformation(), G);
236  }
237 
238  void torus_mesh::adapt(const getfem::mesh &original_mesh){
239  clear();
240  GMM_ASSERT1(original_mesh.dim() == 2, "Adapting torus feature must be a 2d mesh");
241  mesh::copy_from(original_mesh);
242  adapt();
243  }
244 
245  void torus_mesh::adapt(){
246 
247  bgeot::node_tab node_tab_copy(pts);
248  this->pts.clear();
249  for(size_type pt = 0; pt < node_tab_copy.size(); ++pt){
250  node_tab_copy[pt].resize(3);
251  this->pts.add_node(node_tab_copy[pt]);
252  }
253 
254  for(size_type i = 0; i < this->convex_tab.size(); ++i){
255  bgeot::pconvex_structure pstructure
256  = bgeot::torus_structure_descriptor(convex_tab[i].cstruct);
257  convex_tab[i].cstruct = pstructure;
258  }
259 
260  for(size_type i = 0; i < this->gtab.size(); ++i){
261  bgeot::pgeometric_trans pgeom_trans = bgeot::torus_geom_trans_descriptor(gtab[i]);
262  gtab[i] = pgeom_trans;
263  }
264  }
265 }
getfem::mesh_fem::enumerate_dof
virtual void enumerate_dof() const
Renumber the degrees of freedom.
Definition: getfem_mesh_fem.cc:321
getfem::mesh::clear
void clear()
Erase the mesh.
Definition: getfem_mesh.cc:273
getfem::fem_interpolation_context::M
const base_matrix & M() const
non tau-equivalent transformation matrix.
Definition: getfem_fem.cc:43
getfem::torus_fem::real_grad_base_value
void real_grad_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the gradient of all components of the base functions at the current point of the fem_interpolati...
Definition: getfem_torus.cc:100
bgeot_torus.h
Provides mesh of torus.
bgeot::geotrans_interpolation_context::xreal
const base_node & xreal() const
coordinates of the current point, in the real convex.
Definition: bgeot_geometric_trans.cc:299
getfem::torus_fem::hess_base_value
void hess_base_value(const base_node &, base_tensor &) const
Give the value of all hessians (on ref.
Definition: getfem_torus.cc:65
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem::mesh_fem::set_finite_element
void set_finite_element(size_type cv, pfem pf)
Set the finite element method of a convex.
Definition: getfem_mesh_fem.cc:127
getfem::fem_precomp
pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt, dal::pstatic_stored_object dep)
Handles precomputations for FEM.
Definition: getfem_fem.cc:4333
getfem::mesh_fem::fem_of_element
virtual pfem fem_of_element(size_type cv) const
Return the basic fem associated with an element (if no fem is associated, the function will crash!...
Definition: getfem_mesh_fem.h:444
dal::del_stored_object
void del_stored_object(const pstatic_stored_object &o, bool ignore_unstored)
Delete an object and the object which depend on it.
Definition: dal_static_stored_objects.cc:369
getfem::convex_radius_estimate
scalar_type APIDECL convex_radius_estimate(bgeot::pgeometric_trans pgt, const base_matrix &pts)
rough estimate of the radius of the convex using the largest eigenvalue of the jacobian of the geomet...
Definition: getfem_mesh.cc:798
getfem::torus_mesh_fem
Mesh fem object that adapts.
Definition: getfem_torus.h:97
getfem::virtual_fem
Base class for finite element description.
Definition: getfem_fem.h:255
bgeot::geotrans_interpolation_context::xref
const base_node & xref() const
coordinates of the current point, in the reference convex.
Definition: bgeot_geometric_trans.cc:291
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
getfem::virtual_fem::add_node
void add_node(const pdof_description &d, const base_node &pt, const dal::bit_vector &faces)
internal function adding a node to an element for the creation of a finite element method.
Definition: getfem_fem.cc:648
getfem::mesh_fem::convex_index
const dal::bit_vector & convex_index() const
Get the set of convexes where a finite element has been assigned.
Definition: getfem_mesh_fem.h:191
getfem::mesh::copy_from
void copy_from(const mesh &m)
Clone a mesh.
Definition: getfem_mesh.cc:487
getfem::torus_fem
Torus fem, the real grad base value is modified to compute radial grad of F/R.
Definition: getfem_torus.h:52
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem::torus_fem::grad_base_value
void grad_base_value(const base_node &, base_tensor &) const
Give the value of all gradients (on ref.
Definition: getfem_torus.cc:62
getfem::mesh_fem::linked_mesh
const mesh & linked_mesh() const
Return a reference to the underlying mesh.
Definition: getfem_mesh_fem.h:279
getfem::torus_mesh::convex_radius_estimate
virtual scalar_type convex_radius_estimate(size_type ic) const
Return an estimate of the convex largest dimension.
Definition: getfem_torus.cc:228
getfem::torus_fem::base_value
void base_value(const base_node &, base_tensor &) const
Give the value of all components of the base functions at the point x of the reference element.
Definition: getfem_torus.cc:59
bgeot::node_tab
Store a set of points, identifying points that are nearer than a certain very small distance.
Definition: bgeot_node_tab.h:52
getfem::torus_fem::real_hess_base_value
void real_hess_base_value(const fem_interpolation_context &, base_tensor &, bool=true) const
Give the hessian of all components of the base functions at the current point of the fem_interpolatio...
Definition: getfem_torus.cc:143
dal::add_stored_object
void add_stored_object(pstatic_stored_object_key k, pstatic_stored_object o, permanence perm)
Add an object with two optional dependencies.
Definition: dal_static_stored_objects.cc:284
getfem::torus_mesh_fem::enumerate_dof
void enumerate_dof(void) const
Renumber the degrees of freedom.
Definition: getfem_torus.cc:211
bgeot::pconvex_structure
std::shared_ptr< const convex_structure > pconvex_structure
Pointer on a convex structure description.
Definition: bgeot_convex_structure.h:54
getfem::mesh
Describe a mesh (collection of convexes (elements) and points).
Definition: getfem_mesh.h:95
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
getfem::xfem_dof
pdof_description xfem_dof(pdof_description p, size_type ind)
Description of a special dof for Xfem.
Definition: getfem_fem.cc:433
getfem::torus_fem::real_base_value
void real_base_value(const fem_interpolation_context &c, base_tensor &t, bool=true) const
Give the value of all components of the base functions at the current point of the fem_interpolation_...
Definition: getfem_torus.cc:69
getfem_torus.h
Provides mesh and mesh fem of torus.