Functions | |
| void | distort (MeshBase &mesh, const Real factor, const bool perturb_boundary=false) |
| void | redistribute (MeshBase &mesh, const FunctionBase< Real > &mapfunc) |
| void | translate (MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.) |
| void | rotate (MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.) |
| void | scale (MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.) |
| void | all_tri (MeshBase &mesh) |
| void | smooth (MeshBase &, unsigned int, Real) |
| void | flatten (MeshBase &mesh) |
| void | change_boundary_id (MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id) |
| void | change_subdomain_id (MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id) |
Detailed Description
Tools for Mesh modification.
- Date:
- 2004
Function Documentation
| void libMesh::MeshTools::Modification::all_tri | ( | MeshBase & | mesh | ) |
Converts the 2D quadrilateral elements of a Mesh into triangular elements. Note: Only works for 2D elements! 3D elements are ignored. Note: Probably won't do the right thing for meshes which have been refined previously.
Definition at line 717 of file mesh_modification.C.
References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::err, libMesh::MeshBase::get_boundary_info(), libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_assert(), libMesh::BoundaryInfo::n_boundary_ids(), libMesh::MeshBase::n_elem(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::Elem::node(), libMesh::MeshBase::node(), libMesh::Elem::parent(), libMesh::Elem::point(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::QUAD4, libMesh::QUAD8, libMesh::QUAD9, libMesh::remote_elem, libMesh::BoundaryInfo::remove(), libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), libMesh::TRI3, libMesh::TRI6, and libMesh::Elem::type().
{
// The number of elements in the original mesh before any additions
// or deletions.
const dof_id_type n_orig_elem = mesh.n_elem();
// We store pointers to the newly created elements in a vector
// until they are ready to be added to the mesh. This is because
// adding new elements on the fly can cause reallocation and invalidation
// of existing iterators.
std::vector<Elem*> new_elements;
new_elements.reserve (2*n_orig_elem);
// If the original mesh has boundary data, we carry that over
// to the new mesh with triangular elements.
const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_ids() > 0);
// Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
std::vector<Elem*> new_bndry_elements;
std::vector<unsigned short int> new_bndry_sides;
std::vector<boundary_id_type> new_bndry_ids;
// Iterate over the elements, splitting QUADS into
// pairs of conforming triangles.
// FIXME: This algorithm does not work on refined grids!
{
MeshBase::element_iterator el = mesh.elements_begin();
const MeshBase::element_iterator end = mesh.elements_end();
for (; el!=end; ++el)
{
Elem* elem = *el;
const ElemType etype = elem->type();
// all_tri currently only works on coarse meshes
libmesh_assert (!elem->parent());
// We split the quads using the shorter of the two diagonals
// to maintain the best angle properties.
bool edge_swap = false;
// True if we actually split the current element.
bool split_elem = false;
// The two new triangular elements we will split the quad into.
Elem* tri0 = NULL;
Elem* tri1 = NULL;
switch (etype)
{
case QUAD4:
{
split_elem = true;
tri0 = new Tri3;
tri1 = new Tri3;
// Check for possible edge swap
if ((elem->point(0) - elem->point(2)).size() <
(elem->point(1) - elem->point(3)).size())
{
tri0->set_node(0) = elem->get_node(0);
tri0->set_node(1) = elem->get_node(1);
tri0->set_node(2) = elem->get_node(2);
tri1->set_node(0) = elem->get_node(0);
tri1->set_node(1) = elem->get_node(2);
tri1->set_node(2) = elem->get_node(3);
}
else
{
edge_swap=true;
tri0->set_node(0) = elem->get_node(0);
tri0->set_node(1) = elem->get_node(1);
tri0->set_node(2) = elem->get_node(3);
tri1->set_node(0) = elem->get_node(1);
tri1->set_node(1) = elem->get_node(2);
tri1->set_node(2) = elem->get_node(3);
}
break;
}
case QUAD8:
{
split_elem = true;
tri0 = new Tri6;
tri1 = new Tri6;
Node* new_node = mesh.add_point( (mesh.node(elem->node(0)) +
mesh.node(elem->node(1)) +
mesh.node(elem->node(2)) +
mesh.node(elem->node(3)) / 4)
);
// Check for possible edge swap
if ((elem->point(0) - elem->point(2)).size() <
(elem->point(1) - elem->point(3)).size())
{
tri0->set_node(0) = elem->get_node(0);
tri0->set_node(1) = elem->get_node(1);
tri0->set_node(2) = elem->get_node(2);
tri0->set_node(3) = elem->get_node(4);
tri0->set_node(4) = elem->get_node(5);
tri0->set_node(5) = new_node;
tri1->set_node(0) = elem->get_node(0);
tri1->set_node(1) = elem->get_node(2);
tri1->set_node(2) = elem->get_node(3);
tri1->set_node(3) = new_node;
tri1->set_node(4) = elem->get_node(6);
tri1->set_node(5) = elem->get_node(7);
}
else
{
edge_swap=true;
tri0->set_node(0) = elem->get_node(3);
tri0->set_node(1) = elem->get_node(0);
tri0->set_node(2) = elem->get_node(1);
tri0->set_node(3) = elem->get_node(7);
tri0->set_node(4) = elem->get_node(4);
tri0->set_node(5) = new_node;
tri1->set_node(0) = elem->get_node(1);
tri1->set_node(1) = elem->get_node(2);
tri1->set_node(2) = elem->get_node(3);
tri1->set_node(3) = elem->get_node(5);
tri1->set_node(4) = elem->get_node(6);
tri1->set_node(5) = new_node;
}
break;
}
case QUAD9:
{
split_elem = true;
tri0 = new Tri6;
tri1 = new Tri6;
// Check for possible edge swap
if ((elem->point(0) - elem->point(2)).size() <
(elem->point(1) - elem->point(3)).size())
{
tri0->set_node(0) = elem->get_node(0);
tri0->set_node(1) = elem->get_node(1);
tri0->set_node(2) = elem->get_node(2);
tri0->set_node(3) = elem->get_node(4);
tri0->set_node(4) = elem->get_node(5);
tri0->set_node(5) = elem->get_node(8);
tri1->set_node(0) = elem->get_node(0);
tri1->set_node(1) = elem->get_node(2);
tri1->set_node(2) = elem->get_node(3);
tri1->set_node(3) = elem->get_node(8);
tri1->set_node(4) = elem->get_node(6);
tri1->set_node(5) = elem->get_node(7);
}
else
{
edge_swap=true;
tri0->set_node(0) = elem->get_node(0);
tri0->set_node(1) = elem->get_node(1);
tri0->set_node(2) = elem->get_node(3);
tri0->set_node(3) = elem->get_node(4);
tri0->set_node(4) = elem->get_node(8);
tri0->set_node(5) = elem->get_node(7);
tri1->set_node(0) = elem->get_node(1);
tri1->set_node(1) = elem->get_node(2);
tri1->set_node(2) = elem->get_node(3);
tri1->set_node(3) = elem->get_node(5);
tri1->set_node(4) = elem->get_node(6);
tri1->set_node(5) = elem->get_node(8);
}
break;
}
// No need to split elements that are already triangles
case TRI3:
case TRI6:
continue;
// Try to ignore non-2D elements for now
default:
{
libMesh::err << "Warning, encountered non-2D element "
<< Utility::enum_to_string<ElemType>(etype)
<< " in MeshTools::Modification::all_tri(), hope that's OK..."
<< std::endl;
}
} // end switch (etype)
if (split_elem)
{
// Be sure the correct ID's are also set for tri0 and
// tri1.
tri0->processor_id() = elem->processor_id();
tri0->subdomain_id() = elem->subdomain_id();
tri1->processor_id() = elem->processor_id();
tri1->subdomain_id() = elem->subdomain_id();
if (mesh_has_boundary_data)
{
for (unsigned short sn=0; sn<elem->n_sides(); ++sn)
{
const std::vector<boundary_id_type>& bc_ids =
mesh.get_boundary_info().boundary_ids(*el, sn);
for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
{
const boundary_id_type b_id = *id_it;
if (b_id != BoundaryInfo::invalid_id)
{
// Add the boundary ID to the list of new boundary ids
new_bndry_ids.push_back(b_id);
// Convert the boundary side information of the old element to
// boundary side information for the new element.
if (!edge_swap)
{
switch (sn)
{
case 0:
{
// New boundary side is Tri 0, side 0
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(0);
break;
}
case 1:
{
// New boundary side is Tri 0, side 1
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(1);
break;
}
case 2:
{
// New boundary side is Tri 1, side 1
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(1);
break;
}
case 3:
{
// New boundary side is Tri 1, side 2
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(2);
break;
}
default:
libmesh_error_msg("Quad4/8/9 cannot have more than 4 sides.");
}
}
else // edge_swap==true
{
switch (sn)
{
case 0:
{
// New boundary side is Tri 0, side 0
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(0);
break;
}
case 1:
{
// New boundary side is Tri 1, side 0
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(0);
break;
}
case 2:
{
// New boundary side is Tri 1, side 1
new_bndry_elements.push_back(tri1);
new_bndry_sides.push_back(1);
break;
}
case 3:
{
// New boundary side is Tri 0, side 2
new_bndry_elements.push_back(tri0);
new_bndry_sides.push_back(2);
break;
}
default:
libmesh_error_msg("Quad4/8/9 cannot have more than 4 sides.");
}
} // end edge_swap==true
} // end if (b_id != BoundaryInfo::invalid_id)
} // end for loop over boundary IDs
} // end for loop over sides
// Remove the original element from the BoundaryInfo structure.
mesh.get_boundary_info().remove(elem);
} // end if (mesh_has_boundary_data)
// On a distributed mesh, we need to preserve remote_elem
// links, since prepare_for_use can't reconstruct them for
// us.
for (unsigned int sn=0; sn<elem->n_sides(); ++sn)
{
if (elem->neighbor(sn) == remote_elem)
{
// Create a remote_elem link on one of the new
// elements corresponding to the link from the old
// element.
if (!edge_swap)
{
switch (sn)
{
case 0:
{
// New remote side is Tri 0, side 0
tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
break;
}
case 1:
{
// New remote side is Tri 0, side 1
tri0->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
break;
}
case 2:
{
// New remote side is Tri 1, side 1
tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
break;
}
case 3:
{
// New remote side is Tri 1, side 2
tri1->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
break;
}
default:
libmesh_error_msg("Quad4/8/9 cannot have more than 4 sides.");
}
}
else // edge_swap==true
{
switch (sn)
{
case 0:
{
// New remote side is Tri 0, side 0
tri0->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
break;
}
case 1:
{
// New remote side is Tri 1, side 0
tri1->set_neighbor(0, const_cast<RemoteElem*>(remote_elem));
break;
}
case 2:
{
// New remote side is Tri 1, side 1
tri1->set_neighbor(1, const_cast<RemoteElem*>(remote_elem));
break;
}
case 3:
{
// New remote side is Tri 0, side 2
tri0->set_neighbor(2, const_cast<RemoteElem*>(remote_elem));
break;
}
default:
libmesh_error_msg("Quad4/8/9 cannot have more than 4 sides.");
}
} // end edge_swap==true
} // end if (elem->neighbor(sn) == remote_elem)
} // end for loop over sides
// Determine new IDs for the split elements which will be
// the same on all processors, therefore keeping the Mesh
// in sync. Note: we offset the new IDs by n_orig_elem to
// avoid overwriting any of the original IDs, this assumes
// they were contiguously-numbered to begin with...
tri0->set_id( n_orig_elem + 2*elem->id() + 0 );
tri1->set_id( n_orig_elem + 2*elem->id() + 1 );
// Add the newly-created triangles to the temporary vector of new elements.
new_elements.push_back(tri0);
new_elements.push_back(tri1);
// Delete the original element
mesh.delete_elem(elem);
} // end if (split_elem)
} // End for loop over elements
} // end scope
// Now, iterate over the new elements vector, and add them each to
// the Mesh.
{
std::vector<Elem*>::iterator el = new_elements.begin();
const std::vector<Elem*>::iterator end = new_elements.end();
for (; el != end; ++el)
mesh.add_elem(*el);
}
if (mesh_has_boundary_data)
{
// By this time, we should have removed all of the original boundary sides
// - except on a hybrid mesh, where we can't "start from a blank slate"! - RHS
// libmesh_assert_equal_to (mesh.get_boundary_info().n_boundary_conds(), 0);
// Clear the boundary info, to be sure and start from a blank slate.
// mesh.get_boundary_info().clear();
// If the old mesh had boundary data, the new mesh better have some.
libmesh_assert_greater (new_bndry_elements.size(), 0);
// We should also be sure that the lengths of the new boundary data vectors
// are all the same.
libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
// Add the new boundary info to the mesh
for (unsigned int s=0; s<new_bndry_elements.size(); ++s)
mesh.get_boundary_info().add_side(new_bndry_elements[s],
new_bndry_sides[s],
new_bndry_ids[s]);
}
// Prepare the newly created mesh for use.
mesh.prepare_for_use(/*skip_renumber =*/ false);
// Let the new_elements and new_bndry_elements vectors go out of scope.
}
| void libMesh::MeshTools::Modification::change_boundary_id | ( | MeshBase & | mesh, |
| const boundary_id_type | old_id, | ||
| const boundary_id_type | new_id | ||
| ) |
Finds any boundary ids that are currently old_id, changes them to new_id
Definition at line 1490 of file mesh_modification.C.
References libMesh::BoundaryInfo::add_edge(), libMesh::BoundaryInfo::add_node(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::BoundaryInfo::edge_boundary_ids(), libMesh::MeshBase::get_boundary_info(), libMesh::Elem::get_node(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), libMesh::Elem::n_edges(), n_nodes, libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::BoundaryInfo::remove(), libMesh::BoundaryInfo::remove_edge(), and libMesh::BoundaryInfo::remove_side().
{
// Only level-0 elements store BCs. Loop over them.
MeshBase::element_iterator el = mesh.level_elements_begin(0);
const MeshBase::element_iterator end_el = mesh.level_elements_end(0);
for (; el != end_el; ++el)
{
Elem *elem = *el;
unsigned int n_nodes = elem->n_nodes();
for (unsigned int n=0; n != n_nodes; ++n)
{
const std::vector<boundary_id_type>& old_ids =
mesh.get_boundary_info().boundary_ids(elem->get_node(n));
if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
{
std::vector<boundary_id_type> new_ids(old_ids);
std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
mesh.get_boundary_info().remove(elem->get_node(n));
mesh.get_boundary_info().add_node(elem->get_node(n), new_ids);
}
}
unsigned int n_edges = elem->n_edges();
for (unsigned short edge=0; edge != n_edges; ++edge)
{
const std::vector<boundary_id_type>& old_ids =
mesh.get_boundary_info().edge_boundary_ids(elem, edge);
if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
{
std::vector<boundary_id_type> new_ids(old_ids);
std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
mesh.get_boundary_info().remove_edge(elem, edge);
mesh.get_boundary_info().add_edge(elem, edge, new_ids);
}
}
unsigned int n_sides = elem->n_sides();
for (unsigned short s=0; s != n_sides; ++s)
{
const std::vector<boundary_id_type>& old_ids =
mesh.get_boundary_info().boundary_ids(elem, s);
if (std::find(old_ids.begin(), old_ids.end(), old_id) != old_ids.end())
{
std::vector<boundary_id_type> new_ids(old_ids);
std::replace(new_ids.begin(), new_ids.end(), old_id, new_id);
mesh.get_boundary_info().remove_side(elem, s);
mesh.get_boundary_info().add_side(elem, s, new_ids);
}
}
}
}
| void libMesh::MeshTools::Modification::change_subdomain_id | ( | MeshBase & | mesh, |
| const subdomain_id_type | old_id, | ||
| const subdomain_id_type | new_id | ||
| ) |
Finds any subdomain ids that are currently old_id, changes them to new_id
Definition at line 1547 of file mesh_modification.C.
References libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), and libMesh::Elem::subdomain_id().
| void libMesh::MeshTools::Modification::distort | ( | MeshBase & | mesh, |
| const Real | factor, | ||
| const bool | perturb_boundary = false |
||
| ) |
Randomly perturb the nodal locations. This function will move each node factor fraction of its minimum neighboring node separation distance. Nodes on the boundary are not moved by default, however they may be by setting the flag perturb_boundary true.
Definition at line 47 of file mesh_modification.C.
References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::libmesh_assert(), std::max(), libMesh::MeshBase::mesh_dimension(), std::min(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::node_ptr(), libMesh::START_LOG(), and libMesh::TypeVector< T >::unit().
{
libmesh_assert (mesh.n_nodes());
libmesh_assert (mesh.n_elem());
libmesh_assert ((factor >= 0.) && (factor <= 1.));
START_LOG("distort()", "MeshTools::Modification");
// First find nodes on the boundary and flag them
// so that we don't move them
// on_boundary holds false (not on boundary) and true (on boundary)
std::vector<bool> on_boundary (mesh.n_nodes(), false);
if (!perturb_boundary) MeshTools::find_boundary_nodes (mesh, on_boundary);
// Now calculate the minimum distance to
// neighboring nodes for each node.
// hmin holds these distances.
std::vector<float> hmin (mesh.n_nodes(),
std::numeric_limits<float>::max());
MeshBase::element_iterator el = mesh.active_elements_begin();
const MeshBase::element_iterator end = mesh.active_elements_end();
for (; el!=end; ++el)
for (unsigned int n=0; n<(*el)->n_nodes(); n++)
hmin[(*el)->node(n)] = std::min(hmin[(*el)->node(n)],
static_cast<float>((*el)->hmin()));
// Now actually move the nodes
{
const unsigned int seed = 123456;
// seed the random number generator.
// We'll loop from 1 to n_nodes on every processor, even those
// that don't have a particular node, so that the pseudorandom
// numbers will be the same everywhere.
std::srand(seed);
// If the node is on the boundary or
// the node is not used by any element (hmin[n]<1.e20)
// then we should not move it.
// [Note: Testing for (in)equality might be wrong
// (different types, namely float and double)]
for (unsigned int n=0; n<mesh.n_nodes(); n++)
if (!on_boundary[n] && (hmin[n] < 1.e20) )
{
// the direction, random but unit normalized
Point dir( static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
(mesh.mesh_dimension() > 1) ?
static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
: 0.,
((mesh.mesh_dimension() == 3) ?
static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX)
: 0.)
);
dir(0) = (dir(0)-.5)*2.;
if (mesh.mesh_dimension() > 1)
dir(1) = (dir(1)-.5)*2.;
if (mesh.mesh_dimension() == 3)
dir(2) = (dir(2)-.5)*2.;
dir = dir.unit();
Node *node = mesh.node_ptr(n);
if (!node)
continue;
(*node)(0) += dir(0)*factor*hmin[n];
if (mesh.mesh_dimension() > 1)
(*node)(1) += dir(1)*factor*hmin[n];
if (mesh.mesh_dimension() == 3)
(*node)(2) += dir(2)*factor*hmin[n];
}
}
// All done
STOP_LOG("distort()", "MeshTools::Modification");
}
| void libMesh::MeshTools::Modification::flatten | ( | MeshBase & | mesh | ) |
Removes all the refinement tree structure of Mesh, leaving only the highest-level (most-refined) elements. This is useful when you want to write out a uniformly-refined grid to be treated later as an initial mesh. Note that many functions in LibMesh assume a conforming (with no hanging nodes) grid exists at some level, so you probably only want to do this on meshes which have been uniformly refined.
Definition at line 1357 of file mesh_modification.C.
References libMesh::MeshBase::active_elements_begin(), libMesh::MeshBase::active_elements_end(), libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_side(), bc_id, libMesh::BoundaryInfo::boundary_ids(), libMesh::Elem::build(), libMesh::MeshBase::delete_elem(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::MeshBase::get_boundary_info(), libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::BoundaryInfo::invalid_id, libMesh::MeshBase::n_active_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::Elem::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMesh::remote_elem, libMesh::DofObject::set_id(), libMesh::Elem::set_neighbor(), libMesh::Elem::set_node(), libMesh::Elem::subdomain_id(), and libMesh::Elem::type().
{
// Algorithm:
// .) For each active element in the mesh: construct a
// copy which is the same in every way *except* it is
// a level 0 element. Store the pointers to these in
// a separate vector. Save any boundary information as well.
// Delete the active element from the mesh.
// .) Loop over all (remaining) elements in the mesh, delete them.
// .) Add the level-0 copies back to the mesh
// Temporary storage for new element pointers
std::vector<Elem*> new_elements;
// BoundaryInfo Storage for element ids, sides, and BC ids
std::vector<Elem*> saved_boundary_elements;
std::vector<boundary_id_type> saved_bc_ids;
std::vector<unsigned short int> saved_bc_sides;
// Reserve a reasonable amt. of space for each
new_elements.reserve(mesh.n_active_elem());
saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
{
MeshBase::element_iterator it = mesh.active_elements_begin();
const MeshBase::element_iterator end = mesh.active_elements_end();
for (; it != end; ++it)
{
Elem* elem = *it;
// Make a new element of the same type
Elem* copy = Elem::build(elem->type()).release();
// Set node pointers (they still point to nodes in the original mesh)
for(unsigned int n=0; n<elem->n_nodes(); n++)
copy->set_node(n) = elem->get_node(n);
// Copy over ids
copy->processor_id() = elem->processor_id();
copy->subdomain_id() = elem->subdomain_id();
// Retain the original element's ID as well, otherwise ParallelMesh will
// try to create one for you...
copy->set_id( elem->id() );
// This element could have boundary info or ParallelMesh
// remote_elem links as well. We need to save the (elem,
// side, bc_id) triples and those links
for (unsigned short s=0; s<elem->n_sides(); s++)
{
if (elem->neighbor(s) == remote_elem)
copy->set_neighbor(s, const_cast<RemoteElem*>(remote_elem));
const std::vector<boundary_id_type>& bc_ids =
mesh.get_boundary_info().boundary_ids(elem,s);
for (std::vector<boundary_id_type>::const_iterator id_it=bc_ids.begin(); id_it!=bc_ids.end(); ++id_it)
{
const boundary_id_type bc_id = *id_it;
if (bc_id != BoundaryInfo::invalid_id)
{
saved_boundary_elements.push_back(copy);
saved_bc_ids.push_back(bc_id);
saved_bc_sides.push_back(s);
}
}
}
// We're done with this element
mesh.delete_elem(elem);
// But save the copy
new_elements.push_back(copy);
}
// Make sure we saved the same number of boundary conditions
// in each vector.
libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
}
// Loop again, delete any remaining elements
{
MeshBase::element_iterator it = mesh.elements_begin();
const MeshBase::element_iterator end = mesh.elements_end();
for (; it != end; ++it)
mesh.delete_elem( *it );
}
// Add the copied (now level-0) elements back to the mesh
{
for (std::vector<Elem*>::iterator it = new_elements.begin();
it != new_elements.end();
++it)
{
#ifndef NDEBUG
dof_id_type orig_id = (*it)->id();
// ugly mid-statement endif to avoid unused variable warnings
Elem* added_elem =
#endif
mesh.add_elem(*it);
#ifndef NDEBUG
dof_id_type added_id = added_elem->id();
#endif
// If the Elem, as it was re-added to the mesh, now has a
// different ID (this is unlikely, so it's just an assert)
// the boundary information will no longer be correct.
libmesh_assert_equal_to (orig_id, added_id);
}
}
// Finally, also add back the saved boundary information
for (unsigned int e=0; e<saved_boundary_elements.size(); ++e)
mesh.get_boundary_info().add_side(saved_boundary_elements[e],
saved_bc_sides[e],
saved_bc_ids[e]);
// Trim unused and renumber nodes and elements
mesh.prepare_for_use(/*skip_renumber =*/ false);
}
| void libMesh::MeshTools::Modification::redistribute | ( | MeshBase & | mesh, |
| const FunctionBase< Real > & | mapfunc | ||
| ) |
Deterministically perturb the nodal locations. This function will move each node from it's current x/y/z coordinates to a new x/y/z coordinate given by the first LIBMESH_DIM components of the specified function mapfunc
Nodes on the boundary are also moved.
Currently, non-vertex nodes are moved in the same way as vertex nodes, according to (newx,newy,newz) = mapfunc(x,y,z). This behavior is often suboptimal for higher order geometries and may be subject to change in future libMesh versions.
Definition at line 137 of file mesh_modification.C.
References libMesh::FunctionBase< Output >::clone(), end, libMesh::libmesh_assert(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), and libMesh::START_LOG().
{
libmesh_assert (mesh.n_nodes());
libmesh_assert (mesh.n_elem());
START_LOG("redistribute()", "MeshTools::Modification");
DenseVector<Real> output_vec(LIBMESH_DIM);
// FIXME - we should thread this later.
UniquePtr<FunctionBase<Real> > myfunc = mapfunc.clone();
MeshBase::node_iterator it = mesh.nodes_begin();
const MeshBase::node_iterator end = mesh.nodes_end();
for (; it != end; ++it)
{
Node *node = *it;
(*myfunc)(*node, output_vec);
(*node)(0) = output_vec(0);
#if LIBMESH_DIM > 1
(*node)(1) = output_vec(1);
#endif
#if LIBMESH_DIM > 2
(*node)(2) = output_vec(2);
#endif
}
// All done
STOP_LOG("redistribute()", "MeshTools::Modification");
}
| void libMesh::MeshTools::Modification::rotate | ( | MeshBase & | mesh, |
| const Real | phi, | ||
| const Real | theta = 0., |
||
| const Real | psi = 0. |
||
| ) |
Rotates the mesh in the xy plane. The rotation is counter-clock-wise (mathematical definition). The angle is in degrees (360 make a full circle) Rotates the mesh in 3D space. Here the standard Euler angles are adopted (http://mathworld.wolfram.com/EulerAngles.html) The angles are in degrees (360 make a full circle)
Definition at line 211 of file mesh_modification.C.
References libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::pi, libMesh::Real, and libMesh::x.
{
#if LIBMESH_DIM == 3
const Real p = -phi/180.*libMesh::pi;
const Real t = -theta/180.*libMesh::pi;
const Real s = -psi/180.*libMesh::pi;
const Real sp = std::sin(p), cp = std::cos(p);
const Real st = std::sin(t), ct = std::cos(t);
const Real ss = std::sin(s), cs = std::cos(s);
// We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
// (equations 6-14 give the entries of the composite transformation matrix).
// The rotations are performed sequentially about the z, x, and z axes, in that order.
// A positive angle yields a counter-clockwise rotation about the axis in question.
const MeshBase::node_iterator nd_end = mesh.nodes_end();
for (MeshBase::node_iterator nd = mesh.nodes_begin();
nd != nd_end; ++nd)
{
const Point pt = **nd;
const Real x = pt(0);
const Real y = pt(1);
const Real z = pt(2);
**nd = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
(-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
( sp*st)*x + (-cp*st)*y + (ct)*z );
}
#else
libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
#endif
}
| void libMesh::MeshTools::Modification::scale | ( | MeshBase & | mesh, |
| const Real | xs, | ||
| const Real | ys = 0., |
||
| const Real | zs = 0. |
||
| ) |
Scales the mesh. The grid points are scaled in the x direction by xs, in the y direction by ys, etc... If only xs is specified then the scaling is assumed uniform in all directions.
Definition at line 247 of file mesh_modification.C.
References libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::Real, and libMesh::MeshBase::spatial_dimension().
Referenced by libMesh::DenseVector< T >::operator*=(), and libMesh::DenseMatrix< T >::operator*=().
{
const Real x_scale = xs;
Real y_scale = ys;
Real z_scale = zs;
if (ys == 0.)
{
libmesh_assert_equal_to (zs, 0.);
y_scale = z_scale = x_scale;
}
// Scale the x coordinate in all dimensions
const MeshBase::node_iterator nd_end = mesh.nodes_end();
for (MeshBase::node_iterator nd = mesh.nodes_begin();
nd != nd_end; ++nd)
(**nd)(0) *= x_scale;
// Only scale the y coordinate in 2 and 3D
if (mesh.spatial_dimension() < 2)
return;
for (MeshBase::node_iterator nd = mesh.nodes_begin();
nd != nd_end; ++nd)
(**nd)(1) *= y_scale;
// Only scale the z coordinate in 3D
if (mesh.spatial_dimension() < 3)
return;
for (MeshBase::node_iterator nd = mesh.nodes_begin();
nd != nd_end; ++nd)
(**nd)(2) *= z_scale;
}
| void libMesh::MeshTools::Modification::smooth | ( | MeshBase & | mesh, |
| unsigned int | n_iterations, | ||
| Real | power | ||
| ) |
Smooth the mesh with a simple Laplace smoothing algorithm. The mesh is smoothed n_iterations times. If the parameter power is 0, each node is moved to the average postition of the neighboring connected nodes. If power > 0, the node positions are weighted by their distance. The positions of higher order nodes, and nodes living in refined elements, are calculated from the vertex positions of their parent nodes. Only works in 2D.
- Date:
- 2005
This implementation assumes every element "side" has only 2 nodes.
Definition at line 1175 of file mesh_modification.C.
References libMesh::TypeVector< T >::add(), libMesh::TypeVector< T >::add_scaled(), libMesh::Elem::build_side(), libMesh::Elem::child(), libMesh::Elem::embedding_matrix(), end, libMesh::MeshTools::find_boundary_nodes(), libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::MeshBase::level_elements_begin(), libMesh::MeshBase::level_elements_end(), libMesh::MeshBase::mesh_dimension(), libMesh::Elem::n_children(), libMesh::MeshTools::n_levels(), libMesh::Elem::n_neighbors(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_nodes(), libMesh::Elem::n_second_order_adjacent_vertices(), libMesh::Elem::n_vertices(), libMesh::Elem::neighbor(), libMesh::MeshBase::node(), libMesh::Elem::parent(), libMesh::Elem::point(), std::pow(), libMesh::Real, libMesh::Elem::second_order_adjacent_vertex(), side, libMesh::TypeVector< T >::size(), and libMesh::MeshTools::weight().
{
libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
/*
* find the boundary nodes
*/
std::vector<bool> on_boundary;
MeshTools::find_boundary_nodes(mesh, on_boundary);
for (unsigned int iter=0; iter<n_iterations; iter++)
{
/*
* loop over the mesh refinement level
*/
unsigned int n_levels = MeshTools::n_levels(mesh);
for (unsigned int refinement_level=0; refinement_level != n_levels;
refinement_level++)
{
// initialize the storage (have to do it on every level to get empty vectors
std::vector<Point> new_positions;
std::vector<Real> weight;
new_positions.resize(mesh.n_nodes());
weight.resize(mesh.n_nodes());
{
/*
* Loop over the elements to calculate new node positions
*/
MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
for (; el != end; ++el)
{
/*
* Constant handle for the element
*/
const Elem* elem = *el;
/*
* We relax all nodes on level 0 first
* If the element is refined (level > 0), we interpolate the
* parents nodes with help of the embedding matrix
*/
if (refinement_level == 0)
{
for (unsigned int s=0; s<elem->n_neighbors(); s++)
{
/*
* Only operate on sides which are on the
* boundary or for which the current element's
* id is greater than its neighbor's.
* Sides get only built once.
*/
if ((elem->neighbor(s) != NULL) &&
(elem->id() > elem->neighbor(s)->id()) )
{
UniquePtr<Elem> side(elem->build_side(s));
Node* node0 = side->get_node(0);
Node* node1 = side->get_node(1);
Real node_weight = 1.;
// calculate the weight of the nodes
if (power > 0)
{
Point diff = (*node0)-(*node1);
node_weight = std::pow( diff.size(), power );
}
const dof_id_type id0 = node0->id(), id1 = node1->id();
new_positions[id0].add_scaled( *node1, node_weight );
new_positions[id1].add_scaled( *node0, node_weight );
weight[id0] += node_weight;
weight[id1] += node_weight;
}
} // element neighbor loop
}
#ifdef LIBMESH_ENABLE_AMR
else // refinement_level > 0
{
/*
* Find the positions of the hanging nodes of refined elements.
* We do this by calculating their position based on the parent
* (one level less refined) element, and the embedding matrix
*/
const Elem* parent = elem->parent();
/*
* find out which child I am
*/
for (unsigned int c=0; c < parent->n_children(); c++)
{
if (parent->child(c) == elem)
{
/*
*loop over the childs (that is, the current elements) nodes
*/
for (unsigned int nc=0; nc < elem->n_nodes(); nc++)
{
/*
* the new position of the node
*/
Point point;
for (unsigned int n=0; n<parent->n_nodes(); n++)
{
/*
* The value from the embedding matrix
*/
const float em_val = parent->embedding_matrix(c,nc,n);
if (em_val != 0.)
point.add_scaled (parent->point(n), em_val);
}
const dof_id_type id = elem->get_node(nc)->id();
new_positions[id] = point;
weight[id] = 1.;
}
} // if parent->child == elem
} // for parent->n_children
} // if element refinement_level
#endif // #ifdef LIBMESH_ENABLE_AMR
} // element loop
/*
* finally reposition the vertex nodes
*/
for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
if (!on_boundary[nid] && weight[nid] > 0.)
mesh.node(nid) = new_positions[nid]/weight[nid];
}
{
/*
* Now handle the additional second_order nodes by calculating
* their position based on the vertex postitions
* we do a second loop over the level elements
*/
MeshBase::element_iterator el = mesh.level_elements_begin(refinement_level);
const MeshBase::element_iterator end = mesh.level_elements_end(refinement_level);
for (; el != end; ++el)
{
/*
* Constant handle for the element
*/
const Elem* elem = *el;
const unsigned int son_begin = elem->n_vertices();
const unsigned int son_end = elem->n_nodes();
for (unsigned int n=son_begin; n<son_end; n++)
{
const unsigned int n_adjacent_vertices =
elem->n_second_order_adjacent_vertices(n);
Point point;
for (unsigned int v=0; v<n_adjacent_vertices; v++)
point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
const dof_id_type id = elem->get_node(n)->id();
mesh.node(id) = point/n_adjacent_vertices;
}
}
}
} // refinement_level loop
} // end iteration
}
| void libMesh::MeshTools::Modification::translate | ( | MeshBase & | mesh, |
| const Real | xt = 0., |
||
| const Real | yt = 0., |
||
| const Real | zt = 0. |
||
| ) |
Translates the mesh. The grid points are translated in the x direction by xt, in the y direction by yt, etc...
Definition at line 174 of file mesh_modification.C.
References libMesh::MeshBase::nodes_begin(), and libMesh::MeshBase::nodes_end().
generated by
