Hi John, I pasted a somewhat smaller version of my code below. By playing
with the refinement steps limit and the commented out init() and reinit()
lines you should be able to reproduce the two exceptions I got. I hope this
clarifies. Thanks for your time!
//-----------------------------------------------------------------------------------------------------------------------------------------------------
#include <iostream>
#include <vector>
#include <cmath>
#include "libmesh/elem.h"
#include "libmesh/equation_systems.h"
#include "libmesh/linear_implicit_system.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_refinement.h"
#include "libmesh/point.h"
#include "libmesh/system_subset_by_subdomain.h"
#include "libmesh/vtk_io.h"
using namespace libMesh;
using namespace std;
class DummyAssembler : public System::Assembly {
public:
void assemble() {
std::cout << "assembler placeholder" << std::endl;
};
};
enum subdomain_type {
IN_DOMAIN = 1,
SURROGATE = 2,
BETWEEN = 3,
OUTSIDE = 4
};
int main(int argc, char ** argv) {
LibMeshInit init(argc, argv);
Mesh mesh(init.comm());
mesh.read("mesh.msh");
EquationSystems equation_systems(mesh);
equation_systems.add_system<LinearImplicitSystem>("MySystem");
LinearImplicitSystem *system =
&(equation_systems.get_system<LinearImplicitSystem>("MySystem"));
DummyAssembler dummyAssembler;
system->attach_assemble_object(dummyAssembler);
const set<subdomain_id_type> id_list{subdomain_type::IN_DOMAIN,
subdomain_type::SURROGATE};
system->add_variable("u", FIRST, LAGRANGE, &id_list);
{
MeshBase::element_iterator el = mesh.active_local_elements_begin();
const MeshBase::element_iterator elEnd =
mesh.active_local_elements_end();
for (; el != elEnd; ++el) {
(*el)->subdomain_id() = subdomain_type::BETWEEN;
}
}
// equation_systems.init();
Point center(0., 0., 0.);
Real radius = 1.;
size_t refinementSteps = 1;
while (refinementSteps < 2) {
vector<bool> is_included(mesh.max_node_id(), true);
// Checks if node lies inside or outside domain
{
MeshBase::node_iterator node = mesh.active_nodes_begin();
const MeshBase::node_iterator end_node =
mesh.active_nodes_end();
bool inside_circle;
size_t i = 0;
for (; node != end_node; ++node) {
inside_circle = ((Point) **node - center).norm() < radius;
is_included[i] = is_included[i] && inside_circle;
++i;
}
}
// Classifies subdomains ('between' is the one crossed by the
boundary)
{
MeshBase::element_iterator el =
mesh.active_local_subdomain_elements_begin(subdomain_type::BETWEEN);
const MeshBase::element_iterator elEnd =
mesh.active_local_subdomain_elements_end(subdomain_type::BETWEEN);
unsigned count;
for (; el != elEnd; ++el) {
count = 0;
for (size_t i = 0; i < (*el)->n_nodes(); ++i) {
count = count + is_included[(*el)->get_node(i)->id()];
}
if (count == (*el)->n_nodes()) {
(*el)->subdomain_id() = subdomain_type::IN_DOMAIN;
} else if (count == 0) {
(*el)->subdomain_id() = subdomain_type::OUTSIDE;
} else {
(*el)->subdomain_id() = subdomain_type::BETWEEN;
}
}
}
// Refine boundary intersection
{
MeshRefinement meshRefinement(mesh);
MeshBase::element_iterator el =
mesh.active_local_subdomain_elements_begin(subdomain_type::BETWEEN);
const MeshBase::element_iterator elEnd =
mesh.active_local_subdomain_elements_end(subdomain_type::BETWEEN);
for (; el != elEnd; ++el) {
(*el)->set_refinement_flag(Elem::REFINE);
}
meshRefinement.refine_elements();
}
// equation_systems.reinit();
++refinementSteps;
}
// final subdomain evaluation
{
MeshBase::element_iterator el =
mesh.active_local_subdomain_elements_begin(subdomain_type::IN_DOMAIN);
const MeshBase::element_iterator elEnd =
mesh.active_local_subdomain_elements_end(subdomain_type::IN_DOMAIN);
bool surrogate;
for (; el != elEnd; ++el) {
surrogate = false;
for (unsigned i = 0; i < (*el)->n_sides(); ++i) {
surrogate = surrogate ||
(*el)->neighbor_ptr(i)->subdomain_id() == subdomain_type::BETWEEN;
}
if (surrogate) {
(*el)->subdomain_id() = subdomain_type::SURROGATE;
}
}
}
equation_systems.init();
SystemSubsetBySubdomain::SubdomainSelectionByList selection(id_list);
SystemSubsetBySubdomain subset(*system, selection);
system->restrict_solve_to(&subset, SUBSET_ZERO);
VTKIO vtk_io(mesh);
vtk_io.write_equation_systems("out.pvtu", equation_systems);
return 0;
}
//-----------------------------------------------------------------------------------------------------------------------------------------------------
2018-06-11 12:33 GMT-03:00 John Peterson <jwpeter...@gmail.com>:
>
>
> On Sun, Jun 10, 2018 at 10:02 PM, Vinicius C. Reis <
> viniciusr...@coc.ufrj.br> wrote:
>
>> Hi again, it sounds possible that the issue is the hanging nodes. if I
>> only
>> do one refinement step no error occurs.
>
>
> Yes, but presumably there are still hanging nodes even after one level of
> refinement is performed, unless you are referring to one *uniform*
> refinement step...
>
>
>> Adding another step, same
>> exception. Since I am hard flagging the elements to be refined, quite a
>> few
>> hanging nodes appear already in the first refinement step. I tried calling
>> init() before start refining and reinit() at the end of each refinement
>> step, but another error occurred (this one I couldn't check the stack, the
>> debugger didn't manage to stop before finishing the execution).
>>
>
> OK, this is what I was going to suggest (calling reinit() between each
> round of refinement), and it should work. It's possible to you are still
> doing something in your hard flagging step that we are not set up to
> handle, a stack track would be helpful, as would sharing the actual code
> you've written.
>
>
>
>> I am not worried if the refinement is propagated through some of the non
>> selected elements, I was actually expecting this to happen. Is there a way
>> to detect hanging nodes so I can set the neighbor elements to be refined
>> as
>> well?
>
>
> Refining neighbors doesn't get rid of hanging nodes, it just makes a
> hanging node with some other neighbor. You can only eliminate hanging nodes
> by doing uniform refinement...
>
> --
> John
>
------------------------------------------------------------------------------
Check out the vibrant tech community on one of the world's most
engaging tech sites, Slashdot.org! http://sdm.link/slashdot
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
