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