Triangulation<dim> triangulation;
- DoFHandler<dim> dof_handler;
- hp::DoFHandler<dim> hpdof_handler;
- FE_Q<dim> fe;
- hp::FECollection<dim> hpfe;
+ DoFHandler<dim> dof_handler;
+ hp::DoFHandler<dim> hpdof_handler;
+ FE_Q<dim> fe;
+ hp::FECollection<dim> hpfe;
// This is the new variable in the main class. We need an object which holds
// a list of constraints to hold the hanging nodes and the boundary
// conditions.
ConstraintMatrix constraints;
- TimerOutput computing_timer;
+ TimerOutput computing_timer;
SparsityPattern sparsity_pattern;
SparseMatrix<double> system_matrix;
-
+
Vector<double> solution;
Vector<double> system_rhs;
};
template <int dim>
Step6<dim>::Step6 ()
:
- dof_handler (triangulation),
- hpdof_handler (triangulation),
- fe (3),
- computing_timer (std::cout,
- TimerOutput::summary,
- TimerOutput::wall_times)
-
+ dof_handler (triangulation),
+ hpdof_handler (triangulation),
+ fe (3),
+ computing_timer (std::cout,
+ TimerOutput::summary,
+ TimerOutput::wall_times)
+
{
- hpfe.push_back (FE_Q<dim>(3));
+ hpfe.push_back (FE_Q<dim>(3));
}
template <int dim>
void Step6<dim>::setup_system ()
{
- computing_timer.enter_section ("distribute");
+ computing_timer.enter_section ("distribute");
dof_handler.distribute_dofs (fe);
- computing_timer.exit_section ("distribute");
+ computing_timer.exit_section ("distribute");
- computing_timer.enter_section ("distribute_hp");
+ computing_timer.enter_section ("distribute_hp");
hpdof_handler.distribute_dofs (hpfe);
- computing_timer.exit_section ("distribute_hp");
+ computing_timer.exit_section ("distribute_hp");
solution.reinit (dof_handler.n_dofs());
system_rhs.reinit (dof_handler.n_dofs());
// from computations on the previous mesh before the last adaptive
// refinement):
constraints.clear ();
- //computing_timer.enter_section ("hanging");
+ //computing_timer.enter_section ("hanging");
DoFTools::make_hanging_node_constraints (dof_handler,
constraints);
- //computing_timer.exit_section ("hanging");
-/* constraints.clear ();
- computing_timer.enter_section ("hanging_hp");
- DoFTools::make_hanging_node_constraints (dof_handler,
- constraints);
- computing_timer.exit_section ("hanging_hp");
-*/
+ //computing_timer.exit_section ("hanging");
+ /* constraints.clear ();
+ computing_timer.enter_section ("hanging_hp");
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+ computing_timer.exit_section ("hanging_hp");
+ */
// Now we are ready to interpolate the ZeroFunction to our boundary with
// indicator 0 (the whole boundary) and store the resulting constraints in
{
system_matrix = 0;
const QGauss<dim> qformula(fe.degree+1);
-
+
FEValues<dim> fe_values (fe, qformula,
update_values | update_gradients |
update_quadrature_points | update_JxW_values);
FullMatrix<double> cell_matrix;
-
+
Vector<double> cell_rhs;
-
+
std::vector<unsigned int> local_dof_indices;
-
+
const Coefficient<dim> coefficient;
std::vector<double> coefficient_values;
-typename DoFHandler<dim>::active_cell_iterator
+ typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
endc = dof_handler.end();
for (; cell!=endc; ++cell)
{
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
cell_matrix = 0;
cell_rhs.reinit (dofs_per_cell);
cell_rhs = 0;
fe_values.reinit (cell);
coefficient_values.resize(fe_values.n_quadrature_points);
-
- coefficient.value_list (fe_values.get_quadrature_points(),
+
+ coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
// const QGauss<dim> quadrature_formula(3);
hp::QCollection<dim> qformulas;
qformulas.push_back(QGauss<dim>(fe.degree+1));
-
+
hp::FEValues<dim> hp_fe_values (hpfe, qformulas,
- update_values | update_gradients |
- update_quadrature_points | update_JxW_values);
+ update_values | update_gradients |
+ update_quadrature_points | update_JxW_values);
FullMatrix<double> cell_matrix;
-
+
Vector<double> cell_rhs;
-
+
std::vector<unsigned int> local_dof_indices;
-
+
const Coefficient<dim> coefficient;
std::vector<double> coefficient_values;
endc = hpdof_handler.end();
for (; cell!=endc; ++cell)
{
- hp_fe_values.reinit (cell);
- const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
+ hp_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = hp_fe_values.get_present_fe_values ();
const unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell;
- cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
+ cell_matrix.reinit (dofs_per_cell, dofs_per_cell);
cell_matrix = 0;
cell_rhs.reinit (dofs_per_cell);
cell_rhs = 0;
// fe_values.reinit (cell);
coefficient_values.resize(fe_values.n_quadrature_points);
-
- coefficient.value_list (fe_values.get_quadrature_points(),
+
+ coefficient.value_list (fe_values.get_quadrature_points(),
coefficient_values);
for (unsigned int q_point=0; q_point<fe_values.n_quadrature_points; ++q_point)
{
Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
- for (unsigned int i=0;i<triangulation.n_active_cells();++i)
+ for (unsigned int i=0; i<triangulation.n_active_cells(); ++i)
estimated_error_per_cell (i)=i*1.0;
-
- /*
+
+ /*
KellyErrorEstimator<dim>::estimate (dof_handler,
QGauss<dim-1>(3),
typename FunctionMap<dim>::type(),
<< std::endl;
std::cout << "setup" << std::endl;
-
+
// computing_timer.enter_section ("setup");
setup_system ();
// computing_timer.exit_section ("setup");
std::cout << "solve" << std::endl;
- solve ();
+ solve ();
// output_results (cycle);
std::cout << "done" << std::endl;
}
// {Q1,Q2,Q3,Q4}
{
- hp::FECollection<dim> fe_collection;
+ hp::FECollection<dim> fe_collection;
fe_collection.push_back (FE_Q<dim>(1));
fe_collection.push_back (FE_Q<dim>(2));
fe_collection.push_back (FE_Q<dim>(3));
// {Q1xQ1, Q2xQ2, Q3xQ4, Q4xQ3}
{
- hp::FECollection<dim> fe_collection;
+ hp::FECollection<dim> fe_collection;
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(1),1,
FE_Q<dim>(1),1));
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(2),1,
// {Q1xQ1, Q3xQ4, Q4xQ3}
{
- hp::FECollection<dim> fe_collection;
+ hp::FECollection<dim> fe_collection;
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(1),1,
FE_Q<dim>(1),1));
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(3),1,
deallog << fe_collection.find_least_face_dominating_fe(fes) << std::endl;
}
- // {0x0, 0x0, Q1x0, 0xQ1}
+ // {0x0, 0x0, Q1x0, 0xQ1}
{
- hp::FECollection<dim> fe_collection;
+ hp::FECollection<dim> fe_collection;
fe_collection.push_back (FESystem<dim>(FE_Nothing<dim>(),1,
FE_Nothing<dim>(),1));
fe_collection.push_back (FESystem<dim>(FE_Nothing<dim>(),1,
else
deallog << ind << std::endl;
}
-
+
// dominating FE_Nothing
- // {0x0, 0x0, Q1x0, 0xQ1}
+ // {0x0, 0x0, Q1x0, 0xQ1}
{
- hp::FECollection<dim> fe_collection;
+ hp::FECollection<dim> fe_collection;
fe_collection.push_back (FESystem<dim>(FE_Nothing<dim>(1,true),1,
FE_Nothing<dim>(1,true),1));
fe_collection.push_back (FESystem<dim>(FE_Nothing<dim>(1,true),1,
// {Q1xQ1,Q1xQ1,Q2xQ1,Q1,Q2}
{
- hp::FECollection<dim> fe_collection;
- fe_collection.push_back (FESystem<dim>(FE_Q<dim>(1),1,
+ hp::FECollection<dim> fe_collection;
+ fe_collection.push_back (FESystem<dim>(FE_Q<dim>(1),1,
FE_Q<dim>(1),1));
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(1),1,
FE_Q<dim>(1),1));
FE_Q<dim>(2),1));
deallog << fe_collection.find_least_face_dominating_fe(fes) << std::endl;
}
-
+
// {Q4xQ4, Q5xQ5, Q3xQ4, Q4xQ3
{
- hp::FECollection<dim> fe_collection;
+ hp::FECollection<dim> fe_collection;
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(4),1,
FE_Q<dim>(4),1));
fe_collection.push_back (FESystem<dim>(FE_Q<dim>(5),1,
else
deallog << ind << std::endl;
}
-
+
// {Q1,Q2,Q4,Q3}
{
hp::FECollection<dim> fe_collection;
// we put this into a namespace to not conflict with stdlib
namespace Testing
{
-int rand(bool reseed=false, int seed=1) throw()
-{
- static int r[32];
- static int k;
- static bool inited=false;
- if (!inited || reseed)
- {
- //srand treats a seed 0 as 1 for some reason
- r[0]=(seed==0)?1:seed;
-
- for (int i=1; i<31; i++)
- {
- r[i] = (16807LL * r[i-1]) % 2147483647;
- if (r[i] < 0)
- r[i] += 2147483647;
- }
- k=31;
- for (int i=31; i<34; i++)
- {
- r[k%32] = r[(k+32-31)%32];
- k=(k+1)%32;
- }
-
- for (int i=34; i<344; i++)
- {
- r[k%32] = r[(k+32-31)%32] + r[(k+32-3)%32];
- k=(k+1)%32;
- }
- inited=true;
- if (reseed==true)
- return 0;// do not generate new no
- }
+ int rand(bool reseed=false, int seed=1) throw()
+ {
+ static int r[32];
+ static int k;
+ static bool inited=false;
+ if (!inited || reseed)
+ {
+ //srand treats a seed 0 as 1 for some reason
+ r[0]=(seed==0)?1:seed;
- r[k%32] = r[(k+32-31)%32] + r[(k+32-3)%32];
- int ret = r[k%32];
- k=(k+1)%32;
- return (unsigned int)ret >> 1;
-}
+ for (int i=1; i<31; i++)
+ {
+ r[i] = (16807LL * r[i-1]) % 2147483647;
+ if (r[i] < 0)
+ r[i] += 2147483647;
+ }
+ k=31;
+ for (int i=31; i<34; i++)
+ {
+ r[k%32] = r[(k+32-31)%32];
+ k=(k+1)%32;
+ }
+
+ for (int i=34; i<344; i++)
+ {
+ r[k%32] = r[(k+32-31)%32] + r[(k+32-3)%32];
+ k=(k+1)%32;
+ }
+ inited=true;
+ if (reseed==true)
+ return 0;// do not generate new no
+ }
+
+ r[k%32] = r[(k+32-31)%32] + r[(k+32-3)%32];
+ int ret = r[k%32];
+ k=(k+1)%32;
+ return (unsigned int)ret >> 1;
+ }
// reseed our random number generator
-void srand(int seed) throw()
-{
- rand(true, seed);
-}
+ void srand(int seed) throw()
+ {
+ rand(true, seed);
+ }
}
*/
#define check_solver_within_range(SOLVER_COMMAND, CONTROL_COMMAND, MIN_ALLOWED, MAX_ALLOWED) \
-{ \
- const unsigned int previous_depth = deallog.depth_file(0); \
- SOLVER_COMMAND; \
- deallog.depth_file(previous_depth); \
- const unsigned int steps = CONTROL_COMMAND; \
- if (steps >= MIN_ALLOWED && steps <= MAX_ALLOWED) \
- { \
- deallog << "Solver stopped within " << MIN_ALLOWED << " - " \
- << MAX_ALLOWED << " iterations" << std::endl; \
- } \
- else \
- { \
- deallog << "Solver stopped after " << steps << " iterations" \
- << std::endl; \
- } \
-}
+ { \
+ const unsigned int previous_depth = deallog.depth_file(0); \
+ SOLVER_COMMAND; \
+ deallog.depth_file(previous_depth); \
+ const unsigned int steps = CONTROL_COMMAND; \
+ if (steps >= MIN_ALLOWED && steps <= MAX_ALLOWED) \
+ { \
+ deallog << "Solver stopped within " << MIN_ALLOWED << " - " \
+ << MAX_ALLOWED << " iterations" << std::endl; \
+ } \
+ else \
+ { \
+ deallog << "Solver stopped after " << steps << " iterations" \
+ << std::endl; \
+ } \
+ }
// ------------------------------ Functions used in initializing subsystems -------------------
*/
inline unsigned int testing_max_num_threads()
{
- return 5;
+ return 5;
}
struct LimitConcurrency
// stageLog->stageInfo->classLog->classInfo[i].id is always -1, so we look
// it up in stageLog->classLog, make sure it has the same number of entries:
Assert(stageLog->stageInfo->classLog->numClasses == stageLog->classLog->numClasses,
- dealii::ExcInternalError());
+ dealii::ExcInternalError());
bool errors = false;
- for (int i=0;i<stageLog->stageInfo->classLog->numClasses;++i)
+ for (int i=0; i<stageLog->stageInfo->classLog->numClasses; ++i)
{
- if (stageLog->stageInfo->classLog->classInfo[i].destructions !=
- stageLog->stageInfo->classLog->classInfo[i].creations)
- {
- errors = true;
- std::cerr << "ERROR: PETSc objects leaking of type '"
- << stageLog->classLog->classInfo[i].name << "'"
- << " with "
- << stageLog->stageInfo->classLog->classInfo[i].creations
- << " creations and only "
- << stageLog->stageInfo->classLog->classInfo[i].destructions
- << " destructions." << std::endl;
- }
+ if (stageLog->stageInfo->classLog->classInfo[i].destructions !=
+ stageLog->stageInfo->classLog->classInfo[i].creations)
+ {
+ errors = true;
+ std::cerr << "ERROR: PETSc objects leaking of type '"
+ << stageLog->classLog->classInfo[i].name << "'"
+ << " with "
+ << stageLog->stageInfo->classLog->classInfo[i].creations
+ << " creations and only "
+ << stageLog->stageInfo->classLog->classInfo[i].destructions
+ << " destructions." << std::endl;
+ }
}
if (errors)
// enable floating point exceptions
feenableexcept(FE_DIVBYZERO|FE_INVALID);
#endif
- }
+ }
} deal_II_enable_fpe;
projects / dealii.git / commitdiff