#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/fe_q.h>
+#include <deal.II/base/logstream.h>
+
/// to generate random numbers
#include <cstdlib>
+#include <fstream>
using namespace dealii;
-int main(){
- /// Create triangulation, DOF handler and finite element space
+int main()
+{
+ std::ofstream logfile ("anisotropic_crash/output");
+ logfile.precision (3);
+ logfile.setf(std::ios::fixed);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ // Create triangulation
Triangulation<2> tri;
- DoFHandler<2> dh( tri );
- FE_Q<2> fe(1);
GridGenerator::hyper_cube( tri );
tri.refine_global( 2 );
- const unsigned n_cycles = 8;
- for( unsigned i=0; i<n_cycles; ++i ){
- /// Construct the finite element space
- dh.distribute_dofs( fe );
- /// For each vertex find the patch
- const unsigned n_dofs = dh.n_dofs(), n_vert = tri.n_used_vertices();
- std::vector< Point<2> > vertices = tri.get_vertices();
- for( unsigned v=0; v< n_vert; ++v ){
- std::vector< DoFHandler<2>::active_cell_iterator > patch
- = GridTools::find_cells_adjacent_to_vertex( dh, v );
- }
- /// loop over cells and randomly refine
- DoFHandler<2>::active_cell_iterator cell = dh.begin_active(), end = dh.end();
- for( ; cell != end; ++cell ){
- int random = rand()%4;
- /// If random is 0 or 1 we cut x or y, resp.
- if( random < 2 )
- cell->set_refine_flag( RefinementCase<2>::cut_axis( random ) );
- /// If random is 2 we refine isotropically
- if( random == 2 )
- cell->set_refine_flag();
- /// If random is 3 we don't refine
+
+ // now do some sort of random, anisotropic
+ // refinement
+ Triangulation<2>::active_cell_iterator cell = tri.begin_active(), end = tri.end();
+ for( ; cell != end; ++cell )
+ {
+ switch (rand()%4)
+ {
+ /// If a randomly drawn
+ /// number is 0 or 1 we
+ /// cut x or y, resp.
+ case 0:
+ cell->set_refine_flag( RefinementCase<2>::cut_axis(0) );
+ break;
+ case 1:
+ cell->set_refine_flag( RefinementCase<2>::cut_axis(1) );
+ break;
+
+ case 2:
+ /// If the number is 2 we
+ /// refine isotropically
+ cell->set_refine_flag();
+ break;
+
+ default:
+ /// If the number is 3 we don't refine
+ ;
+ }
}
- /// refine the mesh
- tri.prepare_coarsening_and_refinement();
- tri.execute_coarsening_and_refinement();
- }
- dh.clear();
- return 0;
+ /// refine the mesh
+ tri.execute_coarsening_and_refinement();
+
+ /// For each vertex find the patch of cells
+ /// that surrounds it
+ for( unsigned v=0; v<tri.n_used_vertices(); ++v )
+ if (tri.get_used_vertices()[v] == true)
+ GridTools::find_cells_adjacent_to_vertex( tri, v );
}