mesh_smooth.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1097 $
7 //LIC//
8 //LIC// $LastChangedDate: 2015-12-17 11:53:17 +0000 (Thu, 17 Dec 2015) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 
31 #ifndef OOMPH_MESH_SMOOTH_HEADER
32 #define OOMPH_MESH_SMOOTH_HEADER
33 
34 #include<fstream>
35 #include<iostream>
36 
37 #include "../linear_elasticity/elasticity_tensor.h"
38 #include "../constitutive/constitutive_laws.h"
39 #include "../solid/solid_traction_elements.h"
40 
41 
42 
43 namespace oomph
44 {
45 
46 
47 //======================================================================
48 /// Helper namespace
49 //======================================================================
50 namespace Helper_namespace_for_mesh_smoothing
51 {
52 
53  /// Poisson's ratio (for smoothing by linear or nonlinear elasticity)
54  double Nu=0.3;
55 
56  /// Young's modulus (for smoothing by linear or nonlinear elasticity)
57  double E=1.0;
58 
59  /// The elasticity tensor (for smoothing by linear elasticity)
61 
62  /// Create constitutive law (for smoothing by nonlinear elasticity)
64 
65  /// Scale for displacement of quadratic boundary (0.0: simplex; 1.0: quadratic)
66  double Scale=0.1;
67 
68  /// Increment for scale factor for displacement of quadratic boundary
69  double Scale_increment=0.1;
70 
71 }
72 
73 //////////////////////////////////////////////////////////////////
74 //////////////////////////////////////////////////////////////////
75 //////////////////////////////////////////////////////////////////
76 
77 
78 
79 //====================================================================
80 /// Auxiliary Problem to smooth a SolidMesh by adjusting the internal
81 /// nodal positions via the solution of a nonlinear solid mechanics problem.
82 /// The mesh will typically have been created with an unstructured
83 /// mesh generator that uses a low-order (simplex) representation of the
84 /// element geometry; some of the nodes, typically non-vertex nodes on
85 /// the domain's curvilinear boundaries, were then moved to their new
86 /// position to provide a more accurate representation of the geometry.
87 /// This class should be used to deal with elements that may have
88 /// become inverted during the node motion.
89 /// \n
90 /// \b Important \b assumption:
91 /// - It is important that the Lagrangian coordinates of all nodes still
92 /// indicate their original position, i.e. their position before (some of)
93 /// them were moved to their new position. This is because
94 /// we apply the boundary displacements in small increments.
95 /// .
96 /// Template argument specifies type of element. It must be a pure
97 /// solid mechanics element! This shouldn't cause any problems
98 /// since mesh smoothing operations tend to be performed off-line
99 /// so the mesh may as well be built with pure solid elements
100 /// even if it is ultimately to be used with other element types.
101 /// (This restriction could easily be avoided but would
102 /// require double templating and would generally be messy...)
103 //====================================================================
104 template<class ELEMENT>
106 {
107 
108 public:
109 
110  /// \short Functor to update the nodal positions in SolidMesh pointed to by
111  /// orig_mesh_pt in response to the displacement of some of its
112  /// nodes relative to their original position which must still be indicated
113  /// by the nodes' Lagrangian position. copy_of_mesh_pt must be a deep copy
114  /// of orig_mesh_pt, with the same boundary coordinates etc. This mesh
115  /// is used as workspace and can be deleted afterwards.
116  /// The vector controlled_boundary_id contains the ids of
117  /// the mesh boundaries in orig_mesh_pt whose position is supposed
118  /// to remain fixed (while the other nodes are re-positioned to avoid
119  /// the inversion of elements). The final optional argument
120  /// specifies the max. number of increments in which the mesh
121  /// boundary is deformed.
122  void operator()(SolidMesh* orig_mesh_pt, SolidMesh* copy_of_mesh_pt,
123  const Vector<unsigned>& controlled_boundary_id,
124  const unsigned& max_steps=100000000)
125  {
126  // Dummy doc_info
127  DocInfo doc_info;
128  doc_info.disable_doc();
130  copy_of_mesh_pt,
131  controlled_boundary_id,
132  doc_info,
133  max_steps);
134  }
135 
136 
137 
138  /// \short Functor to update the nodal positions in SolidMesh pointed to by
139  /// orig_mesh_pt in response to the displacement of some of its
140  /// nodes relative to their original position which must still be indicated
141  /// by the nodes' Lagrangian position. copy_of_mesh_pt must be a deep copy
142  /// of orig_mesh_pt, with the same boundary coordinates etc. This mesh
143  /// is used as workspace and can be deleted afterwards.
144  /// The vector controlled_boundary_id contains the ids of
145  /// the mesh boundaries in orig_mesh_pt whose position is supposed
146  /// to remain fixed (while the other nodes are re-positioned to avoid
147  /// the inversion of elements). The DocInfo allows allows the output
148  /// of the intermediate meshes. The final optional argument
149  /// specifies the max. number of increments in which the mesh
150  /// boundary is deformed.
151  void operator()(SolidMesh* orig_mesh_pt, SolidMesh* copy_of_mesh_pt,
152  const Vector<unsigned>& controlled_boundary_id,
153  DocInfo doc_info,
154  const unsigned& max_steps=100000000)
155  {
156 
157  // Make original mesh available to everyone...
158  Orig_mesh_pt=orig_mesh_pt;
159  Dummy_mesh_pt=copy_of_mesh_pt;
160 
161  unsigned nnode=orig_mesh_pt->nnode();
162  unsigned nbound=orig_mesh_pt->nboundary();
163  unsigned dim=orig_mesh_pt->node_pt(0)->ndim();
164 
165  // Add to problem's collection of sub-meshes
167 
168  // Backup original nodal positions with boundary nodes snapped
169  // into quadratic position; will soon move these back to
170  // undeformed positon and gently move them back towards
171  // their original position
172  unsigned nnod=Orig_mesh_pt->nnode();
173  Orig_node_pos.resize(nnod);
174  for (unsigned j=0;j<nnod;j++)
175  {
176  Orig_node_pos[j].resize(dim);
177  SolidNode* nod_pt=dynamic_cast<SolidNode*>(Orig_mesh_pt->node_pt(j));
178  for (unsigned i=0;i<dim;i++)
179  {
180  Orig_node_pos[j][i]=nod_pt->x(i);
181  }
182  }
183 
184  // Meshes containing the face elements that represent the
185  // quadratic surface
186  Vector<SolidMesh*> quadratic_surface_mesh_pt(nbound);
187 
188  // GeomObject incarnations
189  Vector<MeshAsGeomObject*> quadratic_surface_geom_obj_pt(nbound);
190 
191 
192  // Create FaceElements on original mesh to define
193  //-------------------------------------------------
194  // the desired boundary shape
195  //---------------------------
196 
197  unsigned n=controlled_boundary_id.size();
198  for (unsigned i=0;i<n;i++)
199  {
200  // Get boundary ID
201  unsigned b=controlled_boundary_id[i];
202 
203  // Create mesh for surface elements
204  quadratic_surface_mesh_pt[b]=new SolidMesh;
205 
206  // How many bulk elements are adjacent to boundary b?
207  unsigned n_element = Orig_mesh_pt->nboundary_element(b);
208 
209  // Loop over the bulk elements adjacent to boundary b
210  for(unsigned e=0;e<n_element;e++)
211  {
212  // Get pointer to the bulk element that is adjacent to boundary b
213  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
215 
216  //What is the index of the face of the element e along boundary b
217  int face_index = Orig_mesh_pt->face_index_at_boundary(b,e);
218 
219  // Create new element
221  new SolidTractionElement<ELEMENT>(bulk_elem_pt,face_index);
222 
223  // Add it to the mesh
224  quadratic_surface_mesh_pt[b]->add_element_pt(el_pt);
225 
226  // Specify boundary number
228  }
229 
230  // Create GeomObject incarnation
231  quadratic_surface_geom_obj_pt[b]=
232  new MeshAsGeomObject(quadratic_surface_mesh_pt[b]);
233  }
234 
235 
236  // Now create Lagrange multiplier elements on dummy mesh
237  //-------------------------------------------------------
238  Vector<SolidMesh*> dummy_lagrange_multiplier_mesh_pt(n);
239  for (unsigned i=0;i<n;i++)
240  {
241  // Get boundary ID
242  unsigned b=controlled_boundary_id[i];
243 
244  // Make new mesh
245  dummy_lagrange_multiplier_mesh_pt[i]=new SolidMesh;
246 
247  // How many bulk elements are adjacent to boundary b?
248  unsigned n_element = Dummy_mesh_pt->nboundary_element(b);
249 
250  // Loop over the bulk fluid elements adjacent to boundary b?
251  for(unsigned e=0;e<n_element;e++)
252  {
253  // Get pointer to the bulk fluid element that is adjacent to boundary b
254  ELEMENT* bulk_elem_pt = dynamic_cast<ELEMENT*>(
256 
257  //Find the index of the face of element e along boundary b
258  int face_index = Dummy_mesh_pt->face_index_at_boundary(b,e);
259 
260  // Create new element
263  bulk_elem_pt,face_index);
264 
265  // Add it to the mesh
266  dummy_lagrange_multiplier_mesh_pt[i]->add_element_pt(el_pt);
267 
268  // Set the GeomObject that defines the boundary shape and set
269  // which bulk boundary we are attached to (needed to extract
270  // the boundary coordinate from the bulk nodes)
272  quadratic_surface_geom_obj_pt[b],b);
273  }
274 
275  // Add sub mesh
276  add_sub_mesh(dummy_lagrange_multiplier_mesh_pt[i]);
277  }
278 
279 
280  // Combine the lot
282 
283  oomph_info << "Number of equations for nonlinear smoothing problem: "
284  << assign_eqn_numbers() << std::endl;
285 
286 
287  // Complete the build of the elements so they are fully functional
288  //----------------------------------------------------------------
289  unsigned n_element = Dummy_mesh_pt->nelement();
290  for(unsigned e=0;e<n_element;e++)
291  {
292  // Upcast from GeneralisedElement to the present element
293  ELEMENT* el_pt =
294  dynamic_cast<ELEMENT*>(Dummy_mesh_pt->element_pt(e));
295 
296  // Set the constitutive law for pseudo-elastic mesh deformation
297  el_pt->constitutive_law_pt() =
299 
300  } // end loop over elements
301 
302 
303  //Output initial configuration
304  doc_solution(doc_info);
305  doc_info.number()++;
306 
307  // Initial scale
310 
311 
312  // Increase scale of deformation until full range is reached
313  //----------------------------------------------------------
314  bool done=false;
315  unsigned count=0;
316  while (!done)
317  {
318 
319  // Increase scale
322 
323  // Backup current nodal positions in dummy mesh
324  backup();
325 
326  // Try it...
327  bool success=true;
328  try
329  {
330  // Avoid overshoot
332  {
334  }
335 
336  // Solve
337  newton_solve();
338  }
340  {
341  success=false;
345 
346  // Reset current nodal positions in dummy mesh
347  reset();
348  }
349 
350  //Output solution
351  if (success)
352  {
353  count++;
354  doc_solution(doc_info);
355  doc_info.number()++;
356  if (Helper_namespace_for_mesh_smoothing::Scale>=1.0) done=true;
357  if (count==max_steps)
358  {
359  oomph_info << "Bailing out after " << count << " steps.\n";
360  done=true;
361  }
362  }
363 
364  }
365 
366  oomph_info << "Done with Helper_namespace_for_mesh_smoothing::Scale="
368 
369  // Loop over nodes in actual mesh and assign new position
370  for (unsigned j=0;j<nnode;j++)
371  {
372  // Get nodes
373  Node* orig_node_pt=orig_mesh_pt->node_pt(j);
374  Node* new_node_pt=Dummy_mesh_pt->node_pt(j);
375 
376  // Assign new position
377  for (unsigned i=0;i<dim;i++)
378  {
379  orig_node_pt->x(i)=new_node_pt->x(i);
380  }
381  }
382 
383  // Now re-assign undeformed position
384  orig_mesh_pt->set_lagrangian_nodal_coordinates();
385 
386  // Cleanup
387  //--------
388  n=controlled_boundary_id.size();
389  for (unsigned i=0;i<n;i++)
390  {
391  // Get boundary ID
392  unsigned b=controlled_boundary_id[i];
393 
394  // Kill meshes and GeomObject representations
395  delete quadratic_surface_mesh_pt[b];
396  delete quadratic_surface_geom_obj_pt[b];
397  delete dummy_lagrange_multiplier_mesh_pt[i];
398  }
399 
400  }
401 
402  /// Destructor (empty)
404 
405 
406  /// \short Update nodal positions in main mesh -- also moves the
407  /// nodes of the FaceElements that impose the new position
409  {
410  oomph_info << "Solving nonlinear smoothing problem for scale "
412  unsigned nnod=Orig_mesh_pt->nnode();
413  for (unsigned j=0;j<nnod;j++)
414  {
415  SolidNode* nod_pt=dynamic_cast<SolidNode*>(Orig_mesh_pt->node_pt(j));
416  unsigned dim=nod_pt->ndim();
417  for (unsigned i=0;i<dim;i++)
418  {
419  nod_pt->x(i)=nod_pt->xi(i)+
421  (Orig_node_pos[j][i]-nod_pt->xi(i));
422  }
423  }
424  }
425 
426 
427  /// \short Backup nodal positions in dummy mesh to allow for reset
428  /// after non-convergence of Newton method
429  void backup()
430  {
431  unsigned nnod=Dummy_mesh_pt->nnode();
432  Backup_node_pos.resize(nnod);
433  for (unsigned j=0;j<nnod;j++)
434  {
435  SolidNode* nod_pt=dynamic_cast<SolidNode*>(Dummy_mesh_pt->node_pt(j));
436  unsigned dim=nod_pt->ndim();
437  Backup_node_pos[j].resize(dim);
438  for (unsigned i=0;i<dim;i++)
439  {
440  Backup_node_pos[j][i]=nod_pt->x(i);
441  }
442  }
443  }
444 
445 
446 
447  /// Reset nodal positions in dummy mesh to allow for restart of
448  /// Newton method with reduced increment in Scale
449  void reset()
450  {
451  unsigned nnod=Dummy_mesh_pt->nnode();
452  for (unsigned j=0;j<nnod;j++)
453  {
454  SolidNode* nod_pt=dynamic_cast<SolidNode*>(Dummy_mesh_pt->node_pt(j));
455  unsigned dim=nod_pt->ndim();
456  for (unsigned i=0;i<dim;i++)
457  {
458  nod_pt->x(i)=Backup_node_pos[j][i];
459  }
460  }
461  }
462 
463  /// Doc the solution
464  void doc_solution(DocInfo& doc_info)
465  {
466 
467  // Bail out
468  if (!doc_info.is_doc_enabled()) return;
469 
470  std::ofstream some_file;
471  std::ostringstream filename;
472 
473  // Number of plot points
474  unsigned npts;
475  npts=5;
476 
477  filename << doc_info.directory() << "/smoothing_soln"
478  << doc_info.number() << ".dat";
479 
480  some_file.open(filename.str().c_str());
481  Dummy_mesh_pt->output(some_file,npts);
482  some_file.close();
483 
484  // Check for inverted elements
485  bool mesh_has_inverted_elements;
486  std::ofstream inverted_fluid_elements;
487  filename.str("");
488  filename << doc_info.directory() << "/inverted_elements_during_smoothing"
489  << doc_info.number() << ".dat";
490  some_file.open(filename.str().c_str());
491  Dummy_mesh_pt->check_inverted_elements(mesh_has_inverted_elements,
492  some_file);
493  some_file.close();
494  oomph_info << "Dummy mesh does ";
495  if (!mesh_has_inverted_elements) oomph_info << "not ";
496  oomph_info << "have inverted elements. \n";
497 
498  }
499 
500  private:
501 
502  /// Original nodal positions
504 
505  /// Backup nodal positions
507 
508  /// Bulk original mesh
510 
511  /// Copy of mesh to work on
513 
514 };
515 
516 
517 //////////////////////////////////////////////////////////////////
518 //////////////////////////////////////////////////////////////////
519 //////////////////////////////////////////////////////////////////
520 
521 
522 
523 //====================================================================
524 /// Auxiliary Problem to smooth a SolidMesh by adjusting the internal
525 /// nodal positions by solving a LINEAR solid mechanics problem for the
526 /// nodal displacements between the specified displacements of certain
527 /// pinned nodes (usually located on boundaries). The template
528 /// parameter specifies the linear elasticity element that must have
529 /// the same shape (geometric element type) as the elements contained
530 /// in the mesh that's to be smoothed. So, e.g. for the ten-noded
531 /// three-dimensional tetrahedral TTaylorHoodElement<3>, it would be
532 /// a TLinearElasticityElement<3,3>, etc.
533 /// \b Important \b assumptions:
534 /// - It is important that the Lagrangian coordinates of all nodes
535 /// still indicate their original position, i.e. their position
536 /// before (some of) them were moved to their new position.
537 /// - It is assumed that in its original state, the mesh does not contain
538 /// any inverted elements.
539 /// .
540 //====================================================================
541 template<class LINEAR_ELASTICITY_ELEMENT>
543 {
544 
545  public:
546 
547  /// \short Constructor: Specify SolidMesh whose nodal positions are to
548  /// be adjusted, and set of nodes in that mesh whose position
549  /// are to remain fixed.
550  void operator()(SolidMesh* orig_mesh_pt,
551  std::set<Node*> pinned_nodes)
552  {
553  // Create new mesh and read out node/element numbers from old one
554  mesh_pt()=new Mesh;
555  unsigned nelem=orig_mesh_pt->nelement();
556  unsigned nnode=orig_mesh_pt->nnode();
557 
558  // Have we already created that node?
559  std::map<Node*,Node*> new_node;
560 
561  // Create new elements
562  for (unsigned e=0;e<nelem;e++)
563  {
564 
565  // Make/add new element
566  LINEAR_ELASTICITY_ELEMENT* el_pt=new LINEAR_ELASTICITY_ELEMENT;
567  mesh_pt()->add_element_pt(el_pt);
568 
569  //Set elasticity tensor
570  el_pt->elasticity_tensor_pt()=
572 
573  // Find corresponding original element
574  SolidFiniteElement* orig_elem_pt=
575  dynamic_cast<SolidFiniteElement*>(orig_mesh_pt->finite_element_pt(e));
576  unsigned nnod=orig_elem_pt->nnode();
577 
578  // Create nodes
579  for (unsigned j=0;j<nnod;j++)
580  {
581  // Does it not exist yet?
582  if (new_node[orig_elem_pt->node_pt(j)]==0)
583  {
584  Node* new_nod_pt=mesh_pt()->finite_element_pt(e)->construct_node(j);
585  new_node[orig_elem_pt->node_pt(j)]=new_nod_pt;
586  mesh_pt()->add_node_pt(new_nod_pt);
587  unsigned dim=new_nod_pt->ndim();
588  for (unsigned i=0;i<dim;i++)
589  {
590  // Set new nodal position to be the old one in the
591  // SolidMesh (assumed to contain no inverted elements)
592  new_nod_pt->x(i)=
593  dynamic_cast<SolidNode*>(orig_elem_pt->node_pt(j))->xi(i);
594  }
595  }
596  // It already exists -- copy across
597  else
598  {
600  new_node[orig_elem_pt->node_pt(j)];
601  }
602  }
603  }
604 
605  // Loop over pinned nodes -- pin their positions and assign updated nodal
606  // positions
607  double scale=1.0;
608  for (std::set<Node*>::iterator it=pinned_nodes.begin();
609  it!=pinned_nodes.end();it++)
610  {
611  unsigned dim=(*it)->ndim();
612  for (unsigned i=0;i<dim;i++)
613  {
614  new_node[*it]->pin(i);
615  new_node[*it]->set_value(i,scale*
616  (dynamic_cast<SolidNode*>(*it)->x(i)-
617  dynamic_cast<SolidNode*>(*it)->xi(i)));
618  }
619  }
620 
621  oomph_info << "Number of equations for smoothing problem: "
622  << assign_eqn_numbers() << std::endl;
623 
624 
625  // Solve
626  newton_solve();
627 
628  // Loop over nodes and assign displacement difference
629  for (unsigned j=0;j<nnode;j++)
630  {
631  // Get nodes
632  SolidNode* orig_node_pt=orig_mesh_pt->node_pt(j);
633  Node* new_node_pt=new_node[orig_node_pt];
634 
635  // Assign displacement difference
636  unsigned dim=new_node_pt->ndim();
637  for (unsigned i=0;i<dim;i++)
638  {
639  orig_node_pt->x(i)=orig_node_pt->xi(i)+new_node_pt->value(i);
640  }
641  }
642 
643  // Now re-assign undeformed position
644  orig_mesh_pt->set_lagrangian_nodal_coordinates();
645 
646  // Clean up -- mesh deletes nodes and elements
647  delete mesh_pt();
648  }
649 
650  /// Destructor (empty)
652 
653 };
654 
655 
656 ////////////////////////////////////////////////////////////////////////
657 ////////////////////////////////////////////////////////////////////////
658 ////////////////////////////////////////////////////////////////////////
659 
660 
661 //====================================================================
662 /// Functor to smooth a SolidMesh by adjusting the internal
663 /// nodal positions by solving a Poisson problem for the
664 /// nodal displacements in the interior. The displacements of the specified
665 /// pinned nodes (usually located on boundaries) remain fixed (their
666 /// displacements are computed from the difference between their
667 /// Lagrangian and Eulerian coordinates). The assumptions is
668 /// that the Lagrangian coordinates in the SolidMesh still reflect
669 /// the original nodal positions before the boundary nodes were
670 /// moved.
671 /// \n
672 /// The template parameter specifies the Poisson element that must have
673 /// the same shape (geometric element type) as the elements contained
674 /// in the mesh that's to be smoothed. So, e.g. for the ten-noded
675 /// three-dimensional tetrahedral TTaylorHoodElement<3>, it would be
676 /// a TPoissonElement<3,3>, etc.
677 //====================================================================
678 template<class POISSON_ELEMENT>
680 {
681 
682 public:
683 
684  /// \short Functor: Specify SolidMesh whose nodal positions are to
685  /// be adjusted, and set of nodes in that mesh whose position
686  /// are to remain fixed.
687  void operator()(SolidMesh* orig_mesh_pt,
688  std::set<Node*> pinned_nodes)
689  {
690  // Create new mesh and read out node/element numbers from old one
691  mesh_pt()=new Mesh;
692  unsigned nelem=orig_mesh_pt->nelement();
693  unsigned nnode=orig_mesh_pt->nnode();
694 
695  // Have we already created that node?
696  std::map<Node*,Node*> new_node;
697 
698  // Create new elements
699  for (unsigned e=0;e<nelem;e++)
700  {
701  mesh_pt()->add_element_pt(new POISSON_ELEMENT);
702 
703  // Find corresponding original element
704  SolidFiniteElement* orig_elem_pt=
705  dynamic_cast<SolidFiniteElement*>(orig_mesh_pt->finite_element_pt(e));
706  unsigned nnod=orig_elem_pt->nnode();
707 
708  // Create nodes
709  for (unsigned j=0;j<nnod;j++)
710  {
711  // Does it not exist yet?
712  if (new_node[orig_elem_pt->node_pt(j)]==0)
713  {
714  Node* new_nod_pt=mesh_pt()->finite_element_pt(e)->construct_node(j);
715  new_node[orig_elem_pt->node_pt(j)]=new_nod_pt;
716  mesh_pt()->add_node_pt(new_nod_pt);
717  unsigned dim=new_nod_pt->ndim();
718  for (unsigned i=0;i<dim;i++)
719  {
720  // Set new nodal position to be the old one in the
721  // SolidMesh (assumed to contain no inverted elements)
722  new_nod_pt->x(i)=
723  dynamic_cast<SolidNode*>(orig_elem_pt->node_pt(j))->xi(i);
724  }
725  }
726  // It already exists -- copy across
727  else
728  {
730  new_node[orig_elem_pt->node_pt(j)];
731  }
732  }
733  }
734 
735 
736  // Loop over pinned nodes
737  for (std::set<Node*>::iterator it=pinned_nodes.begin();
738  it!=pinned_nodes.end();it++)
739  {
740  new_node[*it]->pin(0);
741  }
742 
743  oomph_info << "Number of equations for Poisson displacement smoothing: "
744  << assign_eqn_numbers() << std::endl;
745 
746  // Solve separate displacement problems
747  unsigned dim=orig_mesh_pt->node_pt(0)->ndim();
748  for (unsigned i=0;i<dim;i++)
749  {
750  // Loop over nodes and assign displacement difference
751  for (unsigned j=0;j<nnode;j++)
752  {
753  // Get nodes
754  SolidNode* orig_node_pt=orig_mesh_pt->node_pt(j);
755  Node* new_node_pt=new_node[orig_node_pt];
756 
757  // Assign displacement difference
758  new_node_pt->set_value(0,orig_node_pt->x(i)-orig_node_pt->xi(i));
759  }
760 
761  // Solve
762  newton_solve();
763 
764  // Loop over nodes and assign displacement difference
765  for (unsigned j=0;j<nnode;j++)
766  {
767  // Get nodes
768  SolidNode* orig_node_pt=orig_mesh_pt->node_pt(j);
769  Node* new_node_pt=new_node[orig_node_pt];
770 
771  // Assign displacement difference
772  orig_node_pt->x(i)=orig_node_pt->xi(i)+new_node_pt->value(0);
773  }
774  }
775 
776  // Now re-assign undeformed position
777  orig_mesh_pt->set_lagrangian_nodal_coordinates();
778 
779  // Clean up -- mesh deletes nodes and elements
780  delete mesh_pt();
781 
782  }
783 
784 
785 };
786 
787 
788 /////////////////////////////////////////////////////////////////////////
789 /////////////////////////////////////////////////////////////////////////
790 /////////////////////////////////////////////////////////////////////////
791 
792 }
793 
794 #endif
void operator()(SolidMesh *orig_mesh_pt, SolidMesh *copy_of_mesh_pt, const Vector< unsigned > &controlled_boundary_id, DocInfo doc_info, const unsigned &max_steps=100000000)
Functor to update the nodal positions in SolidMesh pointed to by orig_mesh_pt in response to the disp...
Definition: mesh_smooth.h:151
A class to handle errors in the Newton solver.
Definition: problem.h:2820
void newton_solve()
Use Newton method to solve the problem.
Definition: problem.cc:8741
void operator()(SolidMesh *orig_mesh_pt, std::set< Node * > pinned_nodes)
Constructor: Specify SolidMesh whose nodal positions are to be adjusted, and set of nodes in that mes...
Definition: mesh_smooth.h:550
double value(const unsigned &i) const
Return i-th value (dofs or pinned) at this node either directly or via hanging node representation...
Definition: nodes.cc:2314
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
bool is_doc_enabled() const
Are we documenting?
Information for documentation of results: Directory and file number to enable output in the form RESL...
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
The Problem class.
Definition: problem.h:152
SolidMesh * Dummy_mesh_pt
Copy of mesh to work on.
Definition: mesh_smooth.h:512
void doc_solution(DocInfo &doc_info)
Doc the solution.
Definition: mesh_smooth.h:464
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
OomphInfo oomph_info
e
Definition: cfortran.h:575
void set_boundary_number_in_bulk_mesh(const unsigned &b)
Set function for the boundary number in bulk mesh.
Definition: elements.h:4236
void pin(const unsigned &i)
Pin the i-th stored variable.
Definition: nodes.h:383
unsigned & number()
Number used (e.g.) for labeling output files.
void output(std::ostream &outfile)
Output for all elements.
Definition: mesh.cc:1778
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
void add_node_pt(Node *const &node_pt)
Add a (pointer to a) node to the mesh.
Definition: mesh.h:602
SolidNode * node_pt(const unsigned long &n)
Return a pointer to the n-th global SolidNode.
Definition: mesh.h:2294
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
double & xi(const unsigned &i)
Reference to i-th Lagrangian position.
Definition: nodes.h:1740
void check_inverted_elements(bool &mesh_has_inverted_elements, std::ofstream &inverted_element_file)
Check for inverted elements and report outcome in boolean variable. This visits all elements at their...
Definition: mesh.cc:802
double Scale_increment
Increment for scale factor for displacement of quadratic boundary.
Definition: mesh_smooth.h:69
GeneralisedElement *& element_pt(const unsigned long &e)
Return pointer to element e.
Definition: mesh.h:462
virtual Node * construct_node(const unsigned &n)
Construct the local node n and return a pointer to the newly created node object. ...
Definition: elements.h:2380
void operator()(SolidMesh *orig_mesh_pt, std::set< Node * > pinned_nodes)
Functor: Specify SolidMesh whose nodal positions are to be adjusted, and set of nodes in that mesh wh...
Definition: mesh_smooth.h:687
Vector< Vector< double > > Orig_node_pos
Original nodal positions.
Definition: mesh_smooth.h:503
void set_value(const unsigned &i, const double &value_)
Set the i-th stored data value to specified value. The only reason that we require an explicit set fu...
Definition: nodes.h:267
FiniteElement * boundary_element_pt(const unsigned &b, const unsigned &e) const
Return pointer to e-th finite element on boundary b.
Definition: mesh.h:814
void operator()(SolidMesh *orig_mesh_pt, SolidMesh *copy_of_mesh_pt, const Vector< unsigned > &controlled_boundary_id, const unsigned &max_steps=100000000)
Functor to update the nodal positions in SolidMesh pointed to by orig_mesh_pt in response to the disp...
Definition: mesh_smooth.h:122
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
unsigned nboundary_element(const unsigned &b) const
Return number of finite elements that are adjacent to boundary b.
Definition: mesh.h:852
int face_index_at_boundary(const unsigned &b, const unsigned &e) const
For the e-th finite element on boundary b, return int to indicate the face_index of the face adjacent...
Definition: mesh.h:870
void backup()
Backup nodal positions in dummy mesh to allow for reset after non-convergence of Newton method...
Definition: mesh_smooth.h:429
SolidMesh * Orig_mesh_pt
Bulk original mesh.
Definition: mesh_smooth.h:509
double Scale
Scale for displacement of quadratic boundary (0.0: simplex; 1.0: quadratic)
Definition: mesh_smooth.h:66
unsigned add_sub_mesh(Mesh *const &mesh_pt)
Add a submesh to the problem and return its number, i, by which it can be accessed via mesh_pt(i)...
Definition: problem.h:1293
General SolidMesh class.
Definition: mesh.h:2252
void disable_doc()
Disable documentation.
IsotropicElasticityTensor Isotropic_elasticity_tensor(Nu)
The elasticity tensor (for smoothing by linear elasticity)
std::string directory() const
Output directory.
A Class for nodes that deform elastically (i.e. position is an unknown in the problem). The idea is that the Eulerian positions are stored in a Data object and the Lagrangian coordinates are stored in addition. The pointer that addresses the Eulerian positions is set to the pointer to Value in the Data object. Hence, SolidNode uses knowledge of the internal structure of Data and must be a friend of the Data class. In order to allow a mesh to deform via an elastic-style equation in deforming-domain problems, the positions are stored separately from the values, so that elastic problems may be combined with any other type of problem.
Definition: nodes.h:1568
~LinearElasticitySmoothMesh()
Destructor (empty)
Definition: mesh_smooth.h:651
~NonLinearElasticitySmoothMesh()
Destructor (empty)
Definition: mesh_smooth.h:403
Vector< Vector< double > > Backup_node_pos
Backup nodal positions.
Definition: mesh_smooth.h:506
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void set_boundary_shape_geom_object_pt(GeomObject *boundary_shape_geom_object_pt, const unsigned &boundary_number_in_bulk_mesh)
Set GeomObject that specifies the prescribed boundary displacement; GeomObject is assumed to be param...
void add_element_pt(GeneralisedElement *const &element_pt)
Add a (pointer to) an element to the mesh.
Definition: mesh.h:605
void build_global_mesh()
Build the global mesh by combining the all the submeshes. Note: The nodes boundary information refers...
Definition: problem.cc:1472
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
unsigned nboundary() const
Return number of boundaries.
Definition: mesh.h:806
void actions_before_newton_solve()
Update nodal positions in main mesh – also moves the nodes of the FaceElements that impose the new po...
Definition: mesh_smooth.h:408
double Nu
Poisson's ratio (for smoothing by linear or nonlinear elasticity)
Definition: mesh_smooth.h:54
double E
Young's modulus (for smoothing by linear or nonlinear elasticity)
Definition: mesh_smooth.h:57
ConstitutiveLaw * Constitutive_law_pt
Create constitutive law (for smoothing by nonlinear elasticity)
Definition: mesh_smooth.h:63
void set_lagrangian_nodal_coordinates()
Make the current configuration the undeformed one by setting the nodal Lagrangian coordinates to thei...
Definition: mesh.cc:9196
SolidFiniteElement class.
Definition: elements.h:3320
A general mesh class.
Definition: mesh.h:74