16 #include "adapt_solver.h"
17 #include "solver/linear_solver.h"
18 #include "solver/newton_solver.h"
19 #include "solver/picard_solver.h"
20 #include "ogprojection.h"
22 #include "views/base_view.h"
23 #include "views/order_view.h"
24 #include "views/mesh_view.h"
30 AdaptSolverCriterion::AdaptSolverCriterion()
34 AdaptSolverCriterionErrorThreshold::AdaptSolverCriterionErrorThreshold(
double error_tolerance) : AdaptSolverCriterion(), error_threshold(error_threshold)
38 bool AdaptSolverCriterionErrorThreshold::done(
double error,
unsigned short iteration)
40 return error < this->error_threshold;
43 AdaptSolverCriterionFixed::AdaptSolverCriterionFixed(
unsigned short refinement_levels) : AdaptSolverCriterion(), refinement_levels(refinement_levels)
47 bool AdaptSolverCriterionFixed::done(
double error,
unsigned short iteration)
49 return iteration >= this->refinement_levels;
52 template<
typename Scalar,
typename SolverType>
54 : spaces(initial_spaces), wf(wf), error_calculator(error_calculator), stopping_criterion_single_step(stopping_criterion_single_step), selectors(selectors), stopping_criterion_global(stopping_criterion_global)
59 template<
typename Scalar,
typename SolverType>
61 : wf(wf), stopping_criterion_single_step(stopping_criterion_single_step), error_calculator(error_calculator), stopping_criterion_global(stopping_criterion_global)
63 this->spaces.push_back(initial_space);
64 this->selectors.push_back(selector);
69 template<
typename Scalar,
typename SolverType>
72 this->number_of_equations = this->spaces.size();
73 Helpers::check_length(this->selectors, this->number_of_equations);
74 this->solve_method_running =
false;
75 this->visualization =
false;
76 this->wait_for_keypress =
false;
77 this->prev_mat =
nullptr;
78 this->prev_rhs =
nullptr;
79 this->prev_dirichlet_lift_rhs =
nullptr;
80 this->solver =
new SolverType(wf, spaces);
83 template<
typename Scalar,
typename SolverType>
87 delete this->prev_mat;
90 delete this->prev_rhs;
92 if (this->prev_dirichlet_lift_rhs)
93 delete this->prev_dirichlet_lift_rhs;
98 template<
typename Scalar,
typename SolverType>
101 Loggable::set_verbose_output(to_set);
102 if (this->error_calculator)
103 this->error_calculator->set_verbose_output(to_set);
106 template<
typename Scalar,
typename SolverType>
109 this->visualization = on_off;
110 this->wait_for_keypress = wait_for_keypress_on_off;
113 template<
typename Scalar>
117 static int current_iteration;
118 static std::unordered_set<unsigned int>* current_elements_to_reassemble[H2D_MAX_COMPONENTS];
119 static std::unordered_set<int>* current_DOFs_to_reassemble[H2D_MAX_COMPONENTS];
120 static std::vector<SpaceSharedPtr<Scalar> >* current_ref_spaces;
121 static std::vector<SpaceSharedPtr<Scalar> >* current_prev_ref_spaces;
122 static unsigned char current_number_of_equations;
123 static CSCMatrix<Scalar>* current_prev_mat;
124 static Vector<Scalar>* current_prev_rhs;
125 static Vector<Scalar>* current_prev_dirichlet_lift_rhs;
127 static bool use_Dirichlet;
128 static bool** reusable_DOFs;
129 static bool** reusable_Dirichlet;
132 template<
typename Scalar>
135 template<
typename Scalar>
138 template<
typename Scalar>
141 template<
typename Scalar>
144 template<
typename Scalar>
147 template<
typename Scalar>
150 template<
typename Scalar>
153 template<
typename Scalar>
156 template<
typename Scalar>
159 template<
typename Scalar>
162 template<
typename Scalar>
165 template<
typename Scalar>
170 template<
typename Scalar>
171 static bool compare(Scalar a, Scalar b);
174 bool compare(std::complex<double> a, std::complex<double> b)
176 if (a.real() < b.real() && a.imag() < b.imag())
183 bool compare(
double a,
double b)
188 template<
typename Scalar>
189 void get_states_to_reassemble(Traverse::State**& states,
unsigned int& num_states, SparseMatrix<Scalar>* mat, Vector<Scalar>* rhs, Vector<Scalar>* dirichlet_lift_rhs, Scalar*& coeff_vec)
191 int current_iteration = StateReassemblyHelper<Scalar>::current_iteration;
192 if (StateReassemblyHelper<Scalar>::current_iteration == 1)
194 (*StateReassemblyHelper<Scalar>::reusable_DOFs) =
nullptr;
195 (*StateReassemblyHelper<Scalar>::reusable_Dirichlet) =
nullptr;
201 Hermes::Mixins::TimeMeasurable cpu_time;
204 Space<Scalar>** prev_ref_spaces =
new Space<Scalar>*[StateReassemblyHelper<Scalar>::current_number_of_equations];
205 Space<Scalar>** ref_spaces =
new Space<Scalar>*[StateReassemblyHelper<Scalar>::current_number_of_equations];
207 int total_elements_new_spaces = 0;
209 for (
unsigned short i = 0; i < StateReassemblyHelper<Scalar>::current_number_of_equations; i++)
211 prev_ref_spaces[i] = StateReassemblyHelper<Scalar>::current_prev_ref_spaces->at(i).get();
212 ref_spaces[i] = StateReassemblyHelper<Scalar>::current_ref_spaces->at(i).get();
213 total_elements_new_spaces += ref_spaces[i]->get_mesh()->get_num_active_elements();
217 std::unordered_set<unsigned int> newSpace_elements_to_reassemble[H2D_MAX_COMPONENTS];
219 unsigned int prev_ref_system_size = StateReassemblyHelper<Scalar>::current_prev_mat->get_size();
222 int* DOF_to_DOF_map = malloc_with_check<int>(prev_ref_system_size);
223 for (
unsigned int i = 0; i < prev_ref_system_size; i++)
224 DOF_to_DOF_map[i] = -1;
226 (*StateReassemblyHelper<Scalar>::reusable_DOFs) = calloc_with_check<bool>(ref_system_size);
227 (*StateReassemblyHelper<Scalar>::reusable_Dirichlet) = calloc_with_check<bool>(StateReassemblyHelper<Scalar>::current_number_of_equations);
229 AsmList<Scalar> al, al_prev;
231 std::vector<MeshFunctionSharedPtr<Scalar> > dummy_fns;
232 for (
unsigned short i = 0; i < StateReassemblyHelper<Scalar>::current_number_of_equations; i++)
233 dummy_fns.push_back(
new ZeroSolution<Scalar>(StateReassemblyHelper<Scalar>::current_prev_ref_spaces->at(i)->get_mesh()));
234 for (
unsigned short i = 0; i < StateReassemblyHelper<Scalar>::current_number_of_equations; i++)
235 dummy_fns.push_back(
new ZeroSolution<Scalar>(StateReassemblyHelper<Scalar>::current_ref_spaces->at(i)->get_mesh()));
237 int reusable_DOFs_count = 0;
240 Hermes::Mixins::Loggable::Static::info(
"\t Handling Reusing matrix entries on the new Ref. Space:");
244 for (
unsigned char space_i = 0; space_i < StateReassemblyHelper<Scalar>::current_number_of_equations; space_i++)
246 if (StateReassemblyHelper<Scalar>::current_DOFs_to_reassemble[space_i]->find(-1) == StateReassemblyHelper<Scalar>::current_DOFs_to_reassemble[space_i]->end())
247 (*StateReassemblyHelper<Scalar>::reusable_Dirichlet)[space_i] =
true;
251 Traverse trav(StateReassemblyHelper<Scalar>::current_number_of_equations * 2);
252 unsigned int num_states_local;
253 Traverse::State** states_local = trav.get_states(dummy_fns, num_states_local);
254 for (
unsigned int local_state_i = 0; local_state_i < num_states_local; local_state_i++)
256 Traverse::State* current_local_state = states_local[local_state_i];
257 for (
unsigned char space_i = 0; space_i < StateReassemblyHelper<Scalar>::current_number_of_equations; space_i++)
260 Element* prev_ref_element = current_local_state->e[space_i];
261 Element* ref_element = current_local_state->e[space_i + StateReassemblyHelper<Scalar>::current_number_of_equations];
263 bool element_added =
false;
264 if (StateReassemblyHelper<Scalar>::current_elements_to_reassemble[space_i]->find(prev_ref_element->id) != StateReassemblyHelper<Scalar>::current_elements_to_reassemble[space_i]->end())
266 newSpace_elements_to_reassemble[space_i].insert(ref_element->id);
267 element_added =
true;
270 prev_ref_spaces[space_i]->get_element_assembly_list(prev_ref_element, &al_prev);
271 ref_spaces[space_i]->get_element_assembly_list(ref_element, &al);
273 if (!element_added && al_prev.cnt != al.cnt)
275 newSpace_elements_to_reassemble[space_i].insert(ref_element->id);
276 element_added =
true;
280 unsigned short last_matched_index = 0;
281 bool contains_Dirichlet =
false, anything_changed = al_prev.cnt != al.cnt;
282 for (
unsigned short i_al_prev = 0; i_al_prev < al_prev.cnt; i_al_prev++)
285 if (!element_added && al_prev.idx[i_al_prev] != al.idx[i_al_prev])
287 newSpace_elements_to_reassemble[space_i].insert(ref_element->id);
288 element_added =
true;
290 if (al_prev.dof[i_al_prev] < 0)
292 contains_Dirichlet =
true;
295 for (
unsigned short j_al = last_matched_index; j_al < al.cnt; j_al++)
297 if (al.dof[j_al] < 0)
299 contains_Dirichlet =
true;
302 if (al_prev.idx[i_al_prev] == al.idx[j_al] && compare<Scalar>(-HermesEpsilon, al_prev.coef[i_al_prev] - al.coef[j_al]) && compare<Scalar>(al_prev.coef[i_al_prev] - al.coef[j_al], HermesEpsilon))
304 last_matched_index = j_al + 1;
305 if (StateReassemblyHelper<Scalar>::current_DOFs_to_reassemble[space_i]->find(al_prev.dof[i_al_prev]) == StateReassemblyHelper<Scalar>::current_DOFs_to_reassemble[space_i]->end())
307 if (!(*StateReassemblyHelper<Scalar>::reusable_DOFs)[al.dof[j_al]])
308 reusable_DOFs_count++;
309 (*StateReassemblyHelper<Scalar>::reusable_DOFs)[al.dof[j_al]] =
true;
310 DOF_to_DOF_map[al_prev.dof[i_al_prev]] = al.dof[j_al];
318 assert(DOF_to_DOF_map[al_prev.dof[i_al_prev]] == -1);
319 anything_changed =
true;
323 if (anything_changed && contains_Dirichlet)
324 (*StateReassemblyHelper<Scalar>::reusable_Dirichlet)[space_i] =
false;
329 if (StateReassemblyHelper<Scalar>::use_Dirichlet)
331 for (
unsigned int local_state_i = 0; local_state_i < num_states_local; local_state_i++)
333 Traverse::State* current_local_state = states_local[local_state_i];
334 for (
unsigned char space_i = 0; space_i < StateReassemblyHelper<Scalar>::current_number_of_equations; space_i++)
337 if ((*StateReassemblyHelper<Scalar>::reusable_Dirichlet)[space_i])
340 Element* ref_element = current_local_state->e[space_i + StateReassemblyHelper<Scalar>::current_number_of_equations];
342 if (newSpace_elements_to_reassemble[space_i].find(ref_element->id) != newSpace_elements_to_reassemble[space_i].end())
346 Element* prev_ref_element = current_local_state->e[space_i];
349 prev_ref_spaces[space_i]->get_element_assembly_list(prev_ref_element, &al_prev);
350 ref_spaces[space_i]->get_element_assembly_list(ref_element, &al);
353 for (
unsigned short al_i = 0; al_i < al_prev.cnt; al_i++)
355 if (al_prev.dof[al_i] < 0 || al.dof[al_i] < 0)
357 newSpace_elements_to_reassemble[space_i].insert(ref_element->id);
366 Hermes::Mixins::Loggable::Static::info(
"\t No. of DOFs to reassemble: %i / %i total - %2.0f%%.", ref_system_size - reusable_DOFs_count, ref_system_size, ((
float)(ref_system_size - reusable_DOFs_count) / (
float)ref_system_size) * 100.);
367 Hermes::Mixins::Loggable::Static::info(
"\t Search for elements to reassemble: %4.3f s", cpu_time.last());
372 bool* reassembled_0 = (
bool*)calloc(ref_spaces[0]->get_mesh()->get_max_element_id() + 1,
sizeof(bool));
373 bool* reassembled_1 = (
bool*)calloc(ref_spaces[1]->get_mesh()->get_max_element_id() + 1,
sizeof(bool));
375 Traverse::State** new_states = malloc_with_check<Traverse::State*>(num_states,
true);
376 int new_num_states = 0;
377 for (
unsigned int state_i = 0; state_i < num_states; state_i++)
379 for (
unsigned char space_i = 0; space_i < StateReassemblyHelper<Scalar>::current_number_of_equations; space_i++)
381 if (newSpace_elements_to_reassemble[space_i].find(states[state_i]->e[space_i]->
id) != newSpace_elements_to_reassemble[space_i].end())
383 new_states[new_num_states++] = Traverse::State::clone(states[state_i]);
385 reassembled_0[states[state_i]->e[0]->id] =
true;
386 reassembled_1[states[state_i]->e[1]->id] =
true;
394 Views::ScalarView sc_temp_0(
"Reassembled states",
new Views::WinGeom(600, 10, 600, 600));
395 Views::ScalarView sc_temp_1(
"Reassembled states",
new Views::WinGeom(1230, 10, 600, 600));
396 MeshFunctionSharedPtr<double> reassembled_states_function_0(
new ExactSolutionConstantArray<double, bool>(ref_spaces[0]->get_mesh(), reassembled_0));
397 MeshFunctionSharedPtr<double> reassembled_states_function_1(
new ExactSolutionConstantArray<double, bool>(ref_spaces[1]->get_mesh(), reassembled_1));
398 sc_temp_0.show(reassembled_states_function_0);
399 sc_temp_1.show(reassembled_states_function_1);
401 ::free(reassembled_0);
402 ::free(reassembled_1);
405 for (
unsigned int i = 0; i < num_states; i++)
407 free_with_check(states);
408 new_states = realloc_with_check<Traverse::State*>(new_states, new_num_states);
412 Hermes::Mixins::Loggable::Static::info(
"\t No. of states to reassemble: %i / %i total - %2.0f%%.", new_num_states, num_states, ((
float)new_num_states / (
float)num_states) * 100.);
413 num_states = new_num_states;
414 Hermes::Mixins::Loggable::Static::info(
"\t Picking the new states: %4.3f s", cpu_time.last());
418 Scalar* Ax = StateReassemblyHelper<Scalar>::current_prev_mat->get_Ax();
419 int* Ai = StateReassemblyHelper<Scalar>::current_prev_mat->get_Ai();
420 int* Ap = StateReassemblyHelper<Scalar>::current_prev_mat->get_Ap();
421 Scalar* new_coeff_vec =
nullptr;
423 new_coeff_vec = calloc_with_check<Scalar>(ref_system_size,
true);
424 int total_entries = 0;
425 int used_entries = 0;
426 unsigned short current_row_entries;
427 for (
unsigned int i = 0; i < prev_ref_system_size; i++)
429 current_row_entries = Ap[i + 1] - Ap[i];
430 total_entries += current_row_entries;
431 if (DOF_to_DOF_map[i] != -1)
433 for (
int j = 0; j < current_row_entries; j++)
435 if (DOF_to_DOF_map[Ai[Ap[i] + j]] != -1)
437 mat->add(DOF_to_DOF_map[i], DOF_to_DOF_map[Ai[Ap[i] + j]], Ax[Ap[i] + j]);
442 rhs->add(DOF_to_DOF_map[i], StateReassemblyHelper<Scalar>::current_prev_rhs->
get(i));
444 new_coeff_vec[DOF_to_DOF_map[i]] = coeff_vec[i];
447 if (dirichlet_lift_rhs && StateReassemblyHelper<Scalar>::use_Dirichlet)
449 unsigned int running_count = 0;
450 for (
unsigned short space_i = 0; space_i < StateReassemblyHelper<Scalar>::current_number_of_equations; space_i++)
452 if ((*StateReassemblyHelper<Scalar>::reusable_Dirichlet)[space_i])
454 for (
unsigned int dof_i = running_count; dof_i < running_count + prev_ref_spaces[space_i]->get_num_dofs(); dof_i++)
456 if (DOF_to_DOF_map[dof_i] != -1)
457 dirichlet_lift_rhs->add(DOF_to_DOF_map[dof_i], StateReassemblyHelper<Scalar>::current_prev_dirichlet_lift_rhs->
get(dof_i));
460 running_count += prev_ref_spaces[space_i]->get_num_dofs();
466 free_with_check(coeff_vec,
true);
467 coeff_vec = new_coeff_vec;
470 Hermes::Mixins::Loggable::Static::info(
"\t Reused matrix entries: %i / %i total - %2.0f%%.", used_entries, total_entries, ((
float)used_entries / (
float)total_entries) * 100.);
471 Hermes::Mixins::Loggable::Static::info(
"\t Copying the linear system: %4.3f s", cpu_time.last());
472 free_with_check(DOF_to_DOF_map);
475 template<
typename Scalar,
typename SolverType>
476 void AdaptSolver<Scalar, SolverType>::init_solving(
AdaptivityType adaptivityType)
478 this->adaptivity_step = 1;
479 this->solve_method_running =
true;
482 for (
unsigned char i = 0; i < this->spaces.size(); i++)
484 ref_slns.push_back(MeshFunctionSharedPtr<Scalar>(
new Solution<Scalar>));
485 slns.push_back(MeshFunctionSharedPtr<Scalar>(
new Solution<Scalar>));
486 if (this->visualization)
488 this->scalar_views.push_back(
new Views::ScalarView(
"",
new Views::WinGeom(i * 410, 0, 300, 300)));
489 this->scalar_views.back()->get_linearizer()->set_criterion(Views::LinearizerCriterionFixed(4));
490 this->scalar_views.back()->set_title(
"Reference solution #%i", i);
492 this->order_views.push_back(
new Views::OrderView(
"",
new Views::WinGeom(i * 410, 320, 300, 300)));
493 this->order_views.back()->set_title(
"Reference space #%i - orders", i);
495 this->base_views.push_back(
new Views::BaseView<Scalar>(
"",
new Views::WinGeom(i * 410, 640, 300, 300)));
496 this->base_views.back()->get_linearizer()->set_criterion(Views::LinearizerCriterionFixed(4));
500 this->solver->set_verbose_output(this->get_verbose_output());
502 this->solver->set_matrix_export_format(EXPORT_FORMAT_MATLAB_SIMPLE);
505 this->solver->set_rhs_export_format(EXPORT_FORMAT_MATLAB_SIMPLE);
507 this->adaptivity_internal =
new Adapt<Scalar>(spaces, error_calculator, stopping_criterion_single_step);
508 this->adaptivity_internal->set_verbose_output(this->get_verbose_output());
510 StateReassemblyHelper<Scalar>::reusable_DOFs =
new bool*;
511 StateReassemblyHelper<Scalar>::reusable_Dirichlet =
new bool*;
512 this->solver->dp->set_reusable_DOFs(StateReassemblyHelper<Scalar>::reusable_DOFs, StateReassemblyHelper<Scalar>::reusable_Dirichlet);
513 StateReassemblyHelper<Scalar>::use_Dirichlet = this->solver->dp->add_dirichlet_lift;
514 StateReassemblyHelper<Scalar>::current_number_of_equations = this->number_of_equations;
516 this->ref_spaces.clear();
518 if (adaptivityType == hAdaptivity)
520 for (
unsigned char selector_i = 0; selector_i < this->spaces.size(); selector_i++)
521 hOnlySelectors.push_back(
new RefinementSelectors::HOnlySelector<Scalar>());
523 if (adaptivityType == pAdaptivity)
525 for (
unsigned char selector_i = 0; selector_i < this->spaces.size(); selector_i++)
526 pOnlySelectors.push_back(
new RefinementSelectors::POnlySelector<Scalar>());
530 template<
typename Scalar,
typename SolverType>
531 void AdaptSolver<Scalar, SolverType>::deinit_solving()
533 this->solve_method_running =
false;
535 if (this->visualization)
536 for (
unsigned short i = 0; i < this->number_of_equations; i++)
538 delete this->scalar_views[i];
539 delete this->order_views[i];
540 delete this->base_views[i];
543 delete this->adaptivity_internal;
544 delete StateReassemblyHelper<Scalar>::reusable_DOFs;
545 delete StateReassemblyHelper<Scalar>::reusable_Dirichlet;
547 for (
unsigned char selector_i = 0; selector_i < this->hOnlySelectors.size(); selector_i++)
548 delete hOnlySelectors[selector_i];
549 hOnlySelectors.clear();
550 for (
unsigned char selector_i = 0; selector_i < this->pOnlySelectors.size(); selector_i++)
551 delete pOnlySelectors[selector_i];
552 pOnlySelectors.clear();
555 template<
typename Scalar,
typename SolverType>
559 this->init_solving(adaptivityType);
563 this->info(
"\tAdaptSolver step %d:", this->adaptivity_step);
566 prev_ref_spaces = ref_spaces;
568 for (
unsigned char i = 0; i < number_of_equations; i++)
572 typename Space<Scalar>::ReferenceSpaceCreator u_ref_space_creator(spaces[i], adaptivityType == pAdaptivity ? spaces[i]->get_mesh() : ref_mesh, adaptivityType == hAdaptivity ? 0 : 1);
577 this->solver->set_spaces(ref_spaces);
580 StateReassemblyHelper<Scalar>::current_iteration = this->adaptivity_step;
581 StateReassemblyHelper<Scalar>::current_ref_spaces = &ref_spaces;
582 StateReassemblyHelper<Scalar>::current_prev_ref_spaces = &prev_ref_spaces;
583 for (
unsigned char i = 0; i < number_of_equations; i++)
585 StateReassemblyHelper<Scalar>::current_elements_to_reassemble[i] = &this->elements_to_reassemble[i];
586 StateReassemblyHelper<Scalar>::current_DOFs_to_reassemble[i] = &this->DOFs_to_reassemble[i];
591 this->solver->dp->set_reassembled_states_reuse_linear_system_fn(&get_states_to_reassemble<Scalar>);
595 this->solver->sln_vector =
nullptr;
597 this->solver->solve(new_sln_vector);
598 free_with_check(new_sln_vector,
true);
601 this->solver->solve();
603 this->solver->get_linear_matrix_solver()->free();
606 free_with_check<bool>(*StateReassemblyHelper<Scalar>::reusable_DOFs);
607 free_with_check<bool>(*StateReassemblyHelper<Scalar>::reusable_Dirichlet);
611 delete this->prev_mat;
612 this->prev_mat = (Hermes::Algebra::CSCMatrix<Scalar>*)(this->solver->get_jacobian()->duplicate());
613 StateReassemblyHelper<Scalar>::current_prev_mat = this->prev_mat;
616 delete this->prev_rhs;
617 this->prev_rhs = this->solver->get_residual()->duplicate();
618 if (this->solver->dp->add_dirichlet_lift)
619 this->prev_rhs->subtract_vector(this->solver->dp->dirichlet_lift_rhs);
620 StateReassemblyHelper<Scalar>::current_prev_rhs = this->prev_rhs;
622 if (this->solver->dp->add_dirichlet_lift)
624 if (this->prev_dirichlet_lift_rhs)
625 delete this->prev_dirichlet_lift_rhs;
626 this->prev_dirichlet_lift_rhs = this->solver->dp->dirichlet_lift_rhs->duplicate();
627 StateReassemblyHelper<Scalar>::current_prev_dirichlet_lift_rhs = this->prev_dirichlet_lift_rhs;
633 if (this->visualization)
634 this->visualize(ref_spaces);
637 this->info(
"\tProjecting reference solution on coarse mesh.");
641 if (!this->exact_slns.empty())
643 this->info(
"\tCalculating exact error.");
644 this->error_calculator->calculate_errors(slns, exact_slns,
false);
645 double error = this->error_calculator->get_total_error_squared() * 100;
646 this->info(
"\tExact error: %f.", error);
649 this->info(
"\tCalculating error estimate.");
650 this->error_calculator->calculate_errors(slns, ref_slns,
true);
651 double error = this->error_calculator->get_total_error_squared() * 100;
652 this->info(
"\tThe error estimate: %f.", error);
655 if (this->stopping_criterion_global->done(error, this->adaptivity_step))
657 this->deinit_solving();
662 this->info(
"\tAdapting coarse mesh.");
663 total_elements_prev_spaces = 0;
664 for (
unsigned short i = 0; i < number_of_equations; i++)
665 total_elements_prev_spaces += this->spaces[i]->get_mesh()->get_num_active_elements();
667 if (adaptivityType == hpAdaptivity)
668 this->adaptivity_internal->adapt(this->selectors);
669 else if (adaptivityType == hAdaptivity)
670 this->adaptivity_internal->adapt(hOnlySelectors);
672 this->adaptivity_internal->adapt(pOnlySelectors);
674 this->mark_elements_to_reassemble();
677 this->adaptivity_step++;
681 this->deinit_solving();
684 template<
typename Scalar,
typename SolverType>
687 Hermes::Mixins::Loggable::Static::info(
"\t Marking elements to reassemble on the already used Ref. Space:");
691 for (
unsigned char i = 0; i < this->spaces.size(); i++)
693 this->elements_to_reassemble[i].clear();
694 this->DOFs_to_reassemble[i].clear();
697 int total_elements_prev_ref_spaces = 0;
698 for (
unsigned short i = 0; i < number_of_equations; i++)
699 total_elements_prev_ref_spaces += this->ref_spaces[i]->get_mesh()->get_num_active_elements();
702 int valid_elements_to_refine_count = 0;
703 int ref_elements_changed_count = 0;
704 for (
int i = 0; i < this->adaptivity_internal->elements_to_refine_count; i++)
706 ElementToRefine* element_to_refine = &this->adaptivity_internal->elements_to_refine[i];
707 if (!element_to_refine->valid)
709 valid_elements_to_refine_count++;
711 int component = element_to_refine->comp;
712 Element* e = this->ref_spaces[component]->get_mesh()->get_element_fast(element_to_refine->id);
715 this->elements_to_reassemble[component].insert(e->id);
716 ref_elements_changed_count++;
722 this->elements_to_reassemble[component].insert(e->sons[son_i]->id);
723 ref_elements_changed_count++;
727 Hermes::Mixins::Loggable::Static::info(
"\t No. of coarse mesh elements refined: %i / %i total - %2.0f%%.", valid_elements_to_refine_count, total_elements_prev_spaces, ((
float)valid_elements_to_refine_count / (
float)total_elements_prev_spaces) * 100.);
728 Hermes::Mixins::Loggable::Static::info(
"\t No. of fine mesh elements directly changed: %i / %i total - %2.0f%%.", ref_elements_changed_count, total_elements_prev_ref_spaces, ((
float)ref_elements_changed_count / (
float)total_elements_prev_ref_spaces) * 100.);
730 Hermes::Mixins::Loggable::Static::info(
"\t Identify elements that changed: %4.3f s", this->last());
735 for (
unsigned char space_i = 0; space_i < this->number_of_equations; space_i++)
737 for (std::unordered_set<unsigned int>::iterator it = elements_to_reassemble[space_i].begin(); it != elements_to_reassemble[space_i].end(); ++it)
739 Element* e = this->ref_spaces[space_i]->get_mesh()->get_element_fast(*it);
740 this->ref_spaces[space_i]->get_element_assembly_list(e, &al);
741 for (
int j = 0; j < al.cnt; j++)
742 this->DOFs_to_reassemble[space_i].insert(al.dof[j]);
747 Hermes::Mixins::Loggable::Static::info(
"\t Identify DOFs that changed: %4.3f s", this->last());
752 ref_elements_changed_count = 0;
753 for (
unsigned char space_i = 0; space_i < this->number_of_equations; space_i++)
755 for_all_active_elements_fast(this->ref_spaces[space_i]->get_mesh())
757 this->ref_spaces[space_i]->get_element_assembly_list(e, &al);
758 for (
int j = 0; j < al.cnt; j++)
760 if (this->DOFs_to_reassemble[space_i].find(al.dof[j]) != this->DOFs_to_reassemble[space_i].end())
762 this->elements_to_reassemble[space_i].insert(e->id);
763 ref_elements_changed_count++;
770 Hermes::Mixins::Loggable::Static::info(
"\t No. of all fine mesh elements that need to be reassembled: %i / %i total - %2.0f%%.", ref_elements_changed_count, total_elements_prev_ref_spaces, ((
float)ref_elements_changed_count / (
float)total_elements_prev_ref_spaces) * 100.);
772 Hermes::Mixins::Loggable::Static::info(
"\t Looking for elements containing a changed DOF: %4.3f s", this->last());
775 template<
typename Scalar,
typename SolverType>
776 void AdaptSolver<Scalar, SolverType>::visualize(std::vector<SpaceSharedPtr<Scalar> >& ref_spaces)
778 for (
unsigned short i = 0; i < this->number_of_equations; i++)
780 this->scalar_views[i]->show(this->ref_slns[i]);
781 this->order_views[i]->show(ref_spaces[i]);
782 this->base_views[i]->show(ref_spaces[i]);
785 if (this->wait_for_keypress)
789 template<
typename Scalar,
typename SolverType>
792 if (this->solve_method_running)
793 throw Exceptions::Exception(
"AdaptSolver asked for solutions while it was running.");
797 template<
typename Scalar,
typename SolverType>
800 if (this->solve_method_running)
801 throw Exceptions::Exception(
"AdaptSolver asked for a solution while it was running.");
802 return this->slns[index];
805 template<
typename Scalar,
typename SolverType>
808 if (this->solve_method_running)
809 throw Exceptions::Exception(
"AdaptSolver asked for solutions while it was running.");
810 return this->ref_slns;
813 template<
typename Scalar,
typename SolverType>
816 if (this->solve_method_running)
817 throw Exceptions::Exception(
"AdaptSolver asked for a solution while it was running.");
818 return this->ref_slns[index];
821 template<
typename Scalar,
typename SolverType>
824 return this->stopping_criterion_global;
827 template<
typename Scalar,
typename SolverType>
828 void AdaptSolver<Scalar, SolverType>::set_stopping_criterion_global(AdaptSolverCriterion* stopping_criterion_global)
830 this->stopping_criterion_global = stopping_criterion_global;
833 template<
typename Scalar,
typename SolverType>
836 if (this->solve_method_running)
837 throw Exceptions::Exception(
"AdaptSolver asked to change the exact_slns while it was running.");
838 this->exact_slns = exact_slns;
841 template<
typename Scalar,
typename SolverType>
844 if (this->solve_method_running)
845 throw Exceptions::Exception(
"AdaptSolver asked to change the initial spaces while it was running.");
846 Helpers::check_length(this->spaces, spaces);
847 this->spaces = spaces;
849 this->adaptivity_internal->set_spaces(this->spaces);
850 this->solver->set_spaces(this->spaces);
853 template<
typename Scalar,
typename SolverType>
856 if (this->solve_method_running)
857 throw Exceptions::Exception(
"AdaptSolver asked to change the weak formulation while it was running.");
859 this->solver->set_weak_formulation(wf);
862 template<
typename Scalar,
typename SolverType>
863 void AdaptSolver<Scalar, SolverType>::set_error_calculator(ErrorCalculator<Scalar>* error_calculator)
865 this->error_calculator = error_calculator;
868 template<
typename Scalar,
typename SolverType>
869 void AdaptSolver<Scalar, SolverType>::set_stopping_criterion_single_step(AdaptivityStoppingCriterion<Scalar>* stopping_criterion_single_step)
871 this->stopping_criterion_single_step = stopping_criterion_single_step;
872 this->adaptivity_internal->set_strategy(stopping_criterion_single_step);
875 template<
typename Scalar,
typename SolverType>
876 void AdaptSolver<Scalar, SolverType>::set_selectors(std::vector<RefinementSelectors::Selector<Scalar>*> selectors)
878 Helpers::check_length(this->selectors, selectors);
879 this->selectors = selectors;
882 template<
typename Scalar,
typename SolverType>
883 std::vector<SpaceSharedPtr<Scalar> > AdaptSolver<Scalar, SolverType>::get_initial_spaces()
888 template<
typename Scalar,
typename SolverType>
889 WeakFormSharedPtr<Scalar> AdaptSolver<Scalar, SolverType>::get_wf()
894 template<
typename Scalar,
typename SolverType>
900 template<
typename Scalar,
typename SolverType>
903 return this->error_calculator;
906 template<
typename Scalar,
typename SolverType>
907 AdaptivityStoppingCriterion<Scalar>* AdaptSolver<Scalar, SolverType>::get_stopping_criterion_single_step()
909 return this->stopping_criterion_single_step;
912 template<
typename Scalar,
typename SolverType>
913 std::vector<RefinementSelectors::Selector<Scalar>*> AdaptSolver<Scalar, SolverType>::get_selectors()
915 return this->selectors;
918 template<
typename Scalar,
typename SolverType>
919 bool AdaptSolver<Scalar, SolverType>::isOkay()
const
921 this->adaptivity_internal->check();
922 for (
unsigned short i = 0; i < this->number_of_equations; i++)
923 this->spaces[i]->check();
928 template HERMES_API
class AdaptSolver < double, LinearSolver<double> > ;
929 template HERMES_API
class AdaptSolver < std::complex<double>, LinearSolver<std::complex<double > > > ;
AdaptSolver(std::vector< SpaceSharedPtr< Scalar > > initial_spaces, WeakFormSharedPtr< Scalar > wf, ErrorCalculator< Scalar > *error_calculator, AdaptivityStoppingCriterion< Scalar > *stopping_criterion_single_step, std::vector< RefinementSelectors::Selector< Scalar > * > selectors, AdaptSolverCriterion *stopping_criterion_global)
Constructor.
virtual void set_verbose_output(bool to_set)
See Hermes::Mixins::Loggable.
void solve(AdaptivityType adaptivityType)
The main method - solve.
void switch_visualization(bool on_off, bool wait_for_keypress)
Switch visualization on / off.
SolverType * get_solver()
Getters.
Used to pass the instances of Space around.
RefinementType
Possible refinements of an element.
MeshFunctionSharedPtr< Scalar > get_ref_sln(int index)
Get i-th solution.
int get_num_dofs() const
Returns the number of basis functions contained in the space.
File containing ScalarView class.
static void vector_to_solutions(const Scalar *solution_vector, std::vector< SpaceSharedPtr< Scalar > > spaces, std::vector< MeshFunctionSharedPtr< Scalar > > solutions, std::vector< bool > add_dir_lift=std::vector< bool >(), std::vector< int > start_indices=std::vector< int >())
Passes solution components calculated from solution vector as Solutions.
~AdaptSolver()
Destruct this instance.
std::vector< MeshFunctionSharedPtr< Scalar > > get_slns()
Get the solutions.
virtual SpaceSharedPtr< Scalar > create_ref_space(bool assign_dofs=true)
Methods that user calls to get the reference space pointer (has to be properly casted if necessary)...
void set_initial_spaces(std::vector< SpaceSharedPtr< Scalar > >)
Setters.
void set_exact_solutions(std::vector< MeshFunctionSharedPtr< Scalar > > exact_slns)
Add exact solutions for exact solver calculation.
Evaluation of an error between a (coarse) solution and a reference solution.
Class for creating reference space.
static void wait_for_keypress(const char *text=nullptr)
Waits for keypress. Deprecated.
Represents a finite element space over a domain.
A parent of all refinement selectors. Abstract class.
#define H2D_MAX_ELEMENT_SONS
Macros.
MeshFunctionSharedPtr< Scalar > get_sln(int index)
Get i-th solution.
std::vector< MeshFunctionSharedPtr< Scalar > > get_ref_slns()
Get the solutions.
Class for creating reference mesh.
::xsd::cxx::tree::error< char > error
Error condition.
virtual MeshSharedPtr create_ref_mesh()
static void project_global(SpaceSharedPtr< Scalar > space, MatrixFormVol< Scalar > *custom_projection_jacobian, VectorFormVol< Scalar > *custom_projection_residual, Scalar *target_vec)
This method allows to specify your own OG-projection form.
void init()
Common code for the constructors.