1 #include "postprocessing.h"
2 #include "../space/space.h"
3 #include "../function/solution.h"
6 #include "limit_order.h"
7 #include "discrete_problem/discrete_problem_helpers.h"
13 namespace PostProcessing
15 bool VertexBasedLimiter::wider_bounds_on_boundary =
false;
17 template<
typename Scalar>
18 Limiter<Scalar>::Limiter(SpaceSharedPtr<Scalar> space, Scalar* solution_vector) : component_count(1)
20 spaces.push_back(space);
21 this->init(solution_vector);
24 template<
typename Scalar>
25 Limiter<Scalar>::Limiter(std::vector<SpaceSharedPtr<Scalar> > spaces, Scalar* solution_vector) : spaces(spaces), component_count(spaces.size())
27 this->init(solution_vector);
30 template<
typename Scalar>
31 void Limiter<Scalar>::init(Scalar* solution_vector_)
38 Scalar value = solution_vector_[ndof - 1];
41 memcpy(this->solution_vector, solution_vector_,
sizeof(Scalar)* Space<Scalar>::get_num_dofs(this->spaces));
45 throw Exceptions::Exception(
"Wrong combination of space(s) and solution_vector passed to Limiter().");
50 template<
typename Scalar>
51 Limiter<Scalar>::~Limiter()
53 if (this->solution_vector)
54 delete[] this->solution_vector;
57 template<
typename Scalar>
61 warn_if(this->component_count > 1,
"One solution asked from a Limiter, but multiple solutions exist for limiting.");
64 std::vector<MeshFunctionSharedPtr<Scalar> > solutions;
65 solutions.push_back(solution);
66 this->get_solutions(solutions);
67 return solutions.back();
70 template<
typename Scalar>
73 if (solutions.size() != this->component_count)
74 throw Exceptions::Exception(
"Limiter does not have correct number of spaces, solutions.");
76 this->limited_solutions = solutions;
81 if (this->limited_solutions.empty() || this->limited_solutions.size() != this->component_count)
82 throw Exceptions::Exception(
"Limiter failed for unknown reason.");
85 template<
typename Scalar>
88 return this->changed_element_ids;
91 template<
typename Scalar>
97 template<
typename Scalar>
100 return this->solution_vector;
103 template<
typename Scalar>
106 if (solution_vector_)
111 Scalar value = solution_vector_[ndof - 1];
118 throw Exceptions::Exception(
"Wrong combination of space(s) and solution_vector passed to Limiter().");
123 VertexBasedLimiter::VertexBasedLimiter(
SpaceSharedPtr<double> space,
double* solution_vector,
int maximum_polynomial_order)
124 :
Limiter<double>(space, solution_vector)
126 this->init(maximum_polynomial_order);
129 VertexBasedLimiter::VertexBasedLimiter(std::vector<
SpaceSharedPtr<double> > spaces,
double* solution_vector,
int maximum_polynomial_order)
130 : Limiter<double>(spaces, solution_vector)
132 this->init(maximum_polynomial_order);
135 void VertexBasedLimiter::set_p_coarsening_only()
137 this->p_coarsening_only =
true;
140 void VertexBasedLimiter::init(
int maximum_polynomial_order)
142 this->maximum_polynomial_order = maximum_polynomial_order;
143 this->p_coarsening_only =
false;
146 for (
int i = 0; i < this->component_count; i++)
148 if (this->spaces[i]->get_shapeset()->get_id() != 31)
149 throw Exceptions::Exception(
"VertexBasedLimiter designed for L2ShapesetTaylor. Ignore this exception for unforeseen problems.");
152 vertex_min_values =
nullptr;
153 vertex_max_values =
nullptr;
156 this->mixed_derivatives_count = (maximum_polynomial_order)*(maximum_polynomial_order + 1) / 2;
158 this->print_details =
false;
161 VertexBasedLimiter::~VertexBasedLimiter()
163 deallocate_vertex_values();
166 void VertexBasedLimiter::print_detailed_info(
bool print_details_)
168 this->print_details = print_details_;
171 std::vector<std::pair<int, double> > VertexBasedLimiter::get_correction_factors()
const
173 return this->correction_factors;
176 void VertexBasedLimiter::process()
184 prepare_min_max_vertex_values(
true);
190 std::vector<bool> quadratic_correction_done;
192 if (this->get_verbose_output())
193 std::cout <<
"Quadratic correction" << std::endl;
195 for (
int component = 0; component < this->component_count; component++)
197 MeshSharedPtr mesh = this->spaces[component]->get_mesh();
198 if (this->get_verbose_output() && this->component_count > 1)
199 std::cout <<
"Component: " << component << std::endl;
200 for_all_active_elements(e, mesh)
202 bool second_order =
H2D_GET_H_ORDER(this->spaces[component]->get_element_order(e->id)) >= 2 || H2D_GET_V_ORDER(this->spaces[component]->get_element_order(e->id)) >= 2;
205 quadratic_correction_done.push_back(
false);
209 if (this->get_verbose_output())
210 std::cout <<
"Element: " << e->id << std::endl;
212 quadratic_correction_done.push_back(this->impose_quadratic_correction_factor(e, component));
214 if (this->get_verbose_output())
215 std::cout << std::endl;
224 prepare_min_max_vertex_values(
false);
226 if (this->get_verbose_output())
227 std::cout <<
"Linear correction" << std::endl;
230 for (
int component = 0; component < this->component_count; component++)
232 MeshSharedPtr mesh = this->spaces[component]->get_mesh();
233 if (this->get_verbose_output() && this->component_count > 1)
234 std::cout <<
"Component: " << component << std::endl;
235 for_all_active_elements(e, mesh)
237 if (this->get_verbose_output())
238 std::cout <<
"Element: " << e->id << std::endl;
240 bool second_order =
H2D_GET_H_ORDER(this->spaces[component]->get_element_order(e->id)) >= 2 || H2D_GET_V_ORDER(this->spaces[component]->get_element_order(e->id)) >= 2;
241 if (quadratic_correction_done[running_i++] || !second_order)
242 this->impose_linear_correction_factor(e, component);
244 if (this->get_verbose_output())
245 std::cout << std::endl;
253 void VertexBasedLimiter::impose_linear_correction_factor(Element* e,
int component)
255 double correction_factor = std::numeric_limits<double>::infinity();
257 double centroid_value_multiplied = this->get_centroid_value_multiplied(e, component, 0);
258 if (this->get_verbose_output())
259 std::cout << std::endl <<
"center-value: " << centroid_value_multiplied <<
" (" << 0 <<
") ";
261 Solution<double>* sln =
dynamic_cast<Solution<double>*
>(this->limited_solutions[component].get());
263 for (
int i_vertex = 0; i_vertex < e->get_nvert(); i_vertex++)
265 if (this->get_verbose_output())
266 std::cout << std::endl <<
"vertex: " << i_vertex;
268 Node* vertex = e->vn[i_vertex];
272 double vertex_value = sln->get_ref_value_transformed(e, x, y, 0, 0);
274 if (this->get_verbose_output())
275 std::cout <<
"\tvalue: " << vertex_value;
278 if (std::abs(vertex_value - centroid_value_multiplied) < Hermes::HermesSqrtEpsilon)
281 if (this->get_verbose_output())
282 std::cout <<
"\tcenter_value";
285 if (vertex_value > centroid_value_multiplied)
287 fraction = std::min(1., (this->vertex_max_values[component][vertex->id][0] - centroid_value_multiplied) / (vertex_value - centroid_value_multiplied));
288 if (this->get_verbose_output())
289 std::cout <<
"\tmax_value: " << this->vertex_max_values[component][vertex->id][0];
293 fraction = std::min(1., (this->vertex_min_values[component][vertex->id][0] - centroid_value_multiplied) / (vertex_value - centroid_value_multiplied));
294 if (this->get_verbose_output())
295 std::cout <<
"\tmin_value: " << this->vertex_min_values[component][vertex->id][0];
298 correction_factor = std::min(correction_factor, fraction);
300 if (this->get_verbose_output())
301 std::cout << std::endl <<
"correction_factor " << correction_factor << std::endl;
303 if (correction_factor < (1 - 1e-3))
306 this->spaces[component]->get_element_assembly_list(e, &al);
308 for (
int i_basis_fn = 0; i_basis_fn < al.cnt; i_basis_fn++)
310 int order = this->spaces[component]->get_shapeset()->get_order(al.idx[i_basis_fn], e->get_mode());
313 if (this->p_coarsening_only)
314 this->solution_vector[al.dof[i_basis_fn]] = 0.;
316 this->solution_vector[al.dof[i_basis_fn]] *= correction_factor;
320 this->changed_element_ids.push_back(e->id);
321 this->correction_factors.push_back(std::pair<int, double>(1, correction_factor));
325 bool VertexBasedLimiter::impose_quadratic_correction_factor(Element* e,
int component)
327 if (this->get_verbose_output())
328 std::cout <<
"quadratic: ";
330 double correction_factor = std::numeric_limits<double>::infinity();
332 Solution<double>* sln =
dynamic_cast<Solution<double>*
>(this->limited_solutions[component].get());
334 for (
int i_derivative = 1; i_derivative <= 2; i_derivative++)
336 double centroid_value_multiplied = this->get_centroid_value_multiplied(e, component, i_derivative);
337 if (this->get_verbose_output())
338 std::cout << std::endl <<
"center-value: " << centroid_value_multiplied <<
" (" << i_derivative <<
") ";
340 for (
int i_vertex = 0; i_vertex < e->get_nvert(); i_vertex++)
342 if (this->get_verbose_output())
343 std::cout << std::endl <<
"vertex: " << i_vertex;
345 Node* vertex = e->vn[i_vertex];
349 double vertex_value = sln->get_ref_value_transformed(e, x, y, 0, i_derivative);
351 if (this->get_verbose_output())
352 std::cout <<
"\tvalue: " << vertex_value;
355 if (std::abs(vertex_value - centroid_value_multiplied) < Hermes::HermesSqrtEpsilon)
357 if (this->get_verbose_output())
358 std::cout <<
"\tcenter_value";
363 if (vertex_value > centroid_value_multiplied)
365 fraction = std::min(1., (this->vertex_max_values[component][vertex->id][i_derivative] - centroid_value_multiplied) / (vertex_value - centroid_value_multiplied));
366 if (this->get_verbose_output())
367 std::cout <<
"\tmax_value: " << this->vertex_max_values[component][vertex->id][i_derivative];
371 fraction = std::min(1., (this->vertex_min_values[component][vertex->id][i_derivative] - centroid_value_multiplied) / (vertex_value - centroid_value_multiplied));
372 if (this->get_verbose_output())
373 std::cout <<
"\tmin_value: " << this->vertex_min_values[component][vertex->id][i_derivative];
376 correction_factor = std::min(correction_factor, fraction);
380 if (this->get_verbose_output())
381 std::cout << std::endl <<
"correction_factor " << correction_factor << std::endl;
383 if (correction_factor < (1 - 1e-3))
386 this->spaces[component]->get_element_assembly_list(e, &al);
388 for (
int i_basis_fn = 0; i_basis_fn < al.cnt; i_basis_fn++)
390 int order = this->spaces[component]->get_shapeset()->get_order(al.idx[i_basis_fn], e->get_mode());
393 if (this->p_coarsening_only)
394 this->solution_vector[al.dof[i_basis_fn]] = 0.;
396 this->solution_vector[al.dof[i_basis_fn]] *= correction_factor;
400 this->changed_element_ids.push_back(e->id);
401 this->correction_factors.push_back(std::pair<int, double>(2, correction_factor));
408 void VertexBasedLimiter::prepare_min_max_vertex_values(
bool quadratic)
413 deallocate_vertex_values();
414 allocate_vertex_values();
419 for (
int component = 0; component < this->component_count; component++)
421 MeshSharedPtr mesh = this->spaces[component]->get_mesh();
423 for_all_active_elements(e, mesh)
425 for (
int i_vertex = 0; i_vertex < e->get_nvert(); i_vertex++)
427 Node* vertex = e->vn[i_vertex];
428 for (
int i_derivative = (quadratic ? 1 : 0); i_derivative < (quadratic ? this->mixed_derivatives_count : 1); i_derivative++)
430 double element_centroid_value_multiplied = this->get_centroid_value_multiplied(e, component, i_derivative);
431 this->vertex_min_values[component][vertex->id][i_derivative] = std::min(this->vertex_min_values[component][vertex->id][i_derivative], element_centroid_value_multiplied);
432 this->vertex_max_values[component][vertex->id][i_derivative] = std::max(this->vertex_max_values[component][vertex->id][i_derivative], element_centroid_value_multiplied);
433 if (e->en[i_vertex]->bnd && this->wider_bounds_on_boundary)
435 double element_mid_edge_value_multiplied = this->get_edge_midpoint_value_multiplied(e, component, i_derivative, i_vertex);
437 this->vertex_min_values[component][vertex->id][i_derivative] = std::min(this->vertex_min_values[component][vertex->id][i_derivative], element_mid_edge_value_multiplied);
438 this->vertex_max_values[component][vertex->id][i_derivative] = std::max(this->vertex_max_values[component][vertex->id][i_derivative], element_mid_edge_value_multiplied);
440 Node* next_vertex = e->vn[(i_vertex + 1) % e->get_nvert()];
441 this->vertex_min_values[component][next_vertex->id][i_derivative] = std::min(this->vertex_min_values[component][next_vertex->id][i_derivative], element_mid_edge_value_multiplied);
442 this->vertex_max_values[component][next_vertex->id][i_derivative] = std::max(this->vertex_max_values[component][next_vertex->id][i_derivative], element_mid_edge_value_multiplied);
450 double VertexBasedLimiter::get_centroid_value_multiplied(Element* e,
int component,
int mixed_derivative_index)
452 if (mixed_derivative_index > 5)
454 throw Exceptions::MethodNotImplementedException(
"VertexBasedLimiter::get_centroid_value_multiplied only works for first and second derivatives.");
458 Solution<double>* sln =
dynamic_cast<Solution<double>*
>(this->limited_solutions[component].get());
460 if (e->get_mode() == HERMES_MODE_TRIANGLE)
461 result = sln->get_ref_value_transformed(e,
CENTROID_TRI_X, CENTROID_TRI_Y, 0, mixed_derivative_index);
463 result = sln->get_ref_value_transformed(e,
CENTROID_QUAD_X, CENTROID_QUAD_Y, 0, mixed_derivative_index);
468 double VertexBasedLimiter::get_edge_midpoint_value_multiplied(Element* e,
int component,
int mixed_derivative_index,
int edge)
470 if (mixed_derivative_index > 5)
472 throw Exceptions::MethodNotImplementedException(
"VertexBasedLimiter::get_centroid_value_multiplied only works for first and second derivatives.");
476 Solution<double>* sln =
dynamic_cast<Solution<double>*
>(this->limited_solutions[component].get());
482 if (e->get_mode() == HERMES_MODE_TRIANGLE)
499 result = sln->get_ref_value_transformed(e, x, y, 0, mixed_derivative_index);
523 result = sln->get_ref_value_transformed(e, x, y, 0, mixed_derivative_index);
529 void VertexBasedLimiter::allocate_vertex_values()
531 this->vertex_min_values =
new double**[this->component_count];
532 for (
int i = 0; i < this->component_count; i++)
534 this->vertex_min_values[i] =
new double*[this->spaces[i]->get_mesh()->get_max_node_id()];
536 for (
int j = 0; j < this->spaces[i]->get_mesh()->get_max_node_id(); j++)
538 this->vertex_min_values[i][j] =
new double[this->mixed_derivatives_count];
539 for (
int k = 0; k < this->mixed_derivatives_count; k++)
540 this->vertex_min_values[i][j][k] = std::numeric_limits<double>::infinity();
544 this->vertex_max_values =
new double**[this->component_count];
545 for (
int i = 0; i < this->component_count; i++)
547 this->vertex_max_values[i] =
new double*[this->spaces[i]->get_mesh()->get_max_node_id()];
549 for (
int j = 0; j < this->spaces[i]->get_mesh()->get_max_node_id(); j++)
551 this->vertex_max_values[i][j] =
new double[this->mixed_derivatives_count];
552 for (
int k = 0; k < this->mixed_derivatives_count; k++)
553 this->vertex_max_values[i][j][k] = -std::numeric_limits<double>::infinity();
558 void VertexBasedLimiter::deallocate_vertex_values()
560 if (this->vertex_min_values)
562 for (
int i = 0; i < this->component_count; i++)
564 for (
int j = 0; j < this->spaces[i]->get_mesh()->get_max_node_id(); j++)
566 delete[] this->vertex_min_values[i][j];
567 delete[] this->vertex_max_values[i][j];
570 delete[] this->vertex_min_values[i];
571 delete[] this->vertex_max_values[i];
574 delete[] this->vertex_min_values;
575 delete[] this->vertex_max_values;
579 template<
typename Scalar>
582 source_functions.push_back(source_function);
585 template<
typename Scalar>
590 template<
typename Scalar>
593 std::vector<std::string> markers;
594 markers.push_back(marker);
595 return this->calculate(markers);
601 for (
int i = 0; i < this->number_of_integrals; i++)
604 results[i] += results_local[i];
609 void IntegralCalculator<std::complex<double> >::add_results(std::complex<double> * results_local, std::complex<double> * results)
611 for (
int i = 0; i < this->number_of_integrals; i++)
613 #pragma omp critical (integralAddition)
614 results[i] += results_local[i];
618 template<
typename Scalar>
619 VolumetricIntegralCalculator<Scalar>::VolumetricIntegralCalculator(MeshFunctionSharedPtr<Scalar> source_function,
int number_of_integrals) : IntegralCalculator<Scalar>(source_function, number_of_integrals)
623 template<
typename Scalar>
624 VolumetricIntegralCalculator<Scalar>::VolumetricIntegralCalculator(std::vector<MeshFunctionSharedPtr<Scalar> > source_functions,
int number_of_integrals) : IntegralCalculator<Scalar>(source_functions, number_of_integrals)
628 template<
typename Scalar>
632 this->info(
"User markers");
633 for (
unsigned short i = 0; i < markers.size(); i++)
634 this->info(
"\t%s", markers[i].c_str());
637 bool assemble_everywhere =
false;
639 std::vector<int> internal_markers;
640 for (
unsigned short i = 0; i < markers.size(); i++)
642 if (markers[i] == HERMES_ANY)
644 assemble_everywhere =
true;
650 if (internalMarker.valid)
651 internal_markers.push_back(internalMarker.marker);
655 Scalar* result = (Scalar*)calloc(this->number_of_integrals,
sizeof(Scalar));
657 if (!assemble_everywhere && internal_markers.size() == 0)
661 this->info(
"Internal markers");
662 for (
unsigned short i = 0; i < internal_markers.size(); i++)
663 this->info(
"\t%i", internal_markers[i]);
666 int source_functions_size = this->source_functions.size();
668 Traverse trav(source_functions_size);
669 unsigned int num_states;
672 for (
int i = 0; i < source_functions_size; i++)
673 this->source_functions[i]->set_quad_2d(&g_quad_2d_std);
675 #pragma omp parallel num_threads(this->num_threads_used)
680 int thread_number = omp_get_thread_num();
681 int start = (num_states / this->num_threads_used) * thread_number;
682 int end = (num_states / this->num_threads_used) * (thread_number + 1);
683 if (thread_number == this->num_threads_used - 1)
687 this->info(
"Thread %i, states %i - %i", thread_number, start, end);
689 Hermes::Ord* orders = malloc_with_check<Hermes::Ord>(this->number_of_integrals);
690 Func<Hermes::Ord>** func_ord = malloc_with_check<Func<Hermes::Ord>*>(source_functions_size);
692 MeshFunction<Scalar>** source_functions_cloned = malloc_with_check<MeshFunction<Scalar>*>(source_functions_size);
693 for (
int i = 0; i < source_functions_size; i++)
694 source_functions_cloned[i] = this->source_functions[i]->clone();
695 Func<Scalar>** func = malloc_with_check<Func<Scalar>*>(source_functions_size);
697 Scalar* result_thread_local = calloc_with_check<Scalar>(this->number_of_integrals);
698 Scalar* result_local = malloc_with_check<Scalar>(this->number_of_integrals);
699 double jacobian_x_weights[H2D_MAX_INTEGRATION_POINTS_COUNT];
701 for (
int state_i = start; state_i < end; state_i++)
704 this->info(
"Thread %i, state %i", thread_number, state_i);
708 if (!assemble_everywhere)
710 bool target_marker =
false;
711 for (
unsigned short i = 0; i < internal_markers.size(); i++)
713 if (current_state->rep->
marker == internal_markers[i])
715 target_marker =
true;
723 memset(result_local, 0,
sizeof(Scalar)* this->number_of_integrals);
726 for (
int i = 0; i < source_functions_size; i++)
728 if (current_state->e[i])
729 if (current_state->e[i]->
used)
731 source_functions_cloned[i]->set_active_element(current_state->e[i]);
732 source_functions_cloned[i]->set_transform(current_state->sub_idx[i]);
740 for (
int i = 0; i < this->number_of_integrals; i++)
741 orders[i] = Hermes::Ord(0);
743 for (
int i = 0; i < source_functions_size; i++)
745 if (current_state->e[i])
747 if (current_state->e[i]->
used)
748 func_ord[i] =
new Func<Ord>(source_functions_cloned[i]->get_fn_order());
756 this->order(func_ord, orders);
758 for (
int i = 0; i < this->number_of_integrals; i++)
761 int order_int = order.get_order();
764 for (
int i = 0; i < source_functions_size; i++)
766 if (current_state->e[i])
768 if (current_state->e[i]->
used)
769 func[i] =
init_fn(source_functions_cloned[i], order_int);
771 func[i] = init_zero_fn<Scalar>(current_state->rep->get_mode(), order_int);
774 func[i] = init_zero_fn<Scalar>(current_state->rep->get_mode(), order_int);
780 this->integral(n, jacobian_x_weights, func, &geometry, result_local);
782 for (
unsigned short i = 0; i < source_functions_size; i++)
788 for (
unsigned short i = 0; i < this->number_of_integrals; i++)
789 result_thread_local[i] += result_local[i];
792 this->add_results(result_thread_local, result);
793 free_with_check(result_thread_local);
795 for (
unsigned short i = 0; i < source_functions_size; i++)
796 delete source_functions_cloned[i];
797 free_with_check(source_functions_cloned);
798 free_with_check(orders);
799 free_with_check(func_ord);
800 free_with_check(func);
801 free_with_check(result_local);
805 for (
int i = 0; i < num_states; i++)
807 free_with_check(states);
812 template<
typename Scalar>
817 template<
typename Scalar>
818 SurfaceIntegralCalculator<Scalar>::SurfaceIntegralCalculator(std::vector<
MeshFunctionSharedPtr<Scalar> > source_functions,
int number_of_integrals) : IntegralCalculator<Scalar>(source_functions, number_of_integrals)
822 template<
typename Scalar>
826 this->info(
"User markers");
827 for (
unsigned short i = 0; i < markers.size(); i++)
828 this->info(
"\t%s", markers[i].c_str());
831 bool assemble_everywhere =
false;
833 std::vector<int> internal_markers;
834 for (
unsigned short i = 0; i < markers.size(); i++)
836 if (markers[i] == HERMES_ANY)
838 assemble_everywhere =
true;
844 if (internalMarker.valid)
845 internal_markers.push_back(internalMarker.marker);
849 Scalar* result = (Scalar*)calloc(this->number_of_integrals,
sizeof(Scalar));
851 if (!assemble_everywhere && internal_markers.size() == 0)
855 this->info(
"Internal markers");
856 for (
unsigned short i = 0; i < internal_markers.size(); i++)
857 this->info(
"\t%i", internal_markers[i]);
860 int source_functions_size = this->source_functions.size();
861 Traverse trav(source_functions_size);
862 unsigned int num_states;
865 for (
int i = 0; i < source_functions_size; i++)
866 this->source_functions[i]->set_quad_2d(&g_quad_2d_std);
868 #pragma omp parallel num_threads(this->num_threads_used)
873 int thread_number = omp_get_thread_num();
874 int start = (num_states / this->num_threads_used) * thread_number;
875 int end = (num_states / this->num_threads_used) * (thread_number + 1);
876 if (thread_number == this->num_threads_used - 1)
880 this->info(
"Thread %i, states %i - %i", thread_number, start, end);
883 Hermes::Ord* orders = malloc_with_check<Hermes::Ord>(this->number_of_integrals);
884 Func<Hermes::Ord>** func_ord = malloc_with_check<Func<Hermes::Ord>*>(source_functions_size);
886 MeshFunction<Scalar>** source_functions_cloned = malloc_with_check<MeshFunction<Scalar>*>(source_functions_size);
887 for (
int i = 0; i < source_functions_size; i++)
888 source_functions_cloned[i] = this->source_functions[i]->clone();
889 Func<Scalar>** func = malloc_with_check<Func<Scalar>*>(source_functions_size);
891 Scalar* result_thread_local = calloc_with_check<Scalar>(this->number_of_integrals);
892 Scalar* result_local = malloc_with_check<Scalar>(this->number_of_integrals);
893 double jacobian_x_weights[H2D_MAX_INTEGRATION_POINTS_COUNT];
895 for (
int state_i = start; state_i < end; state_i++)
898 this->info(
"Thread %i, state %i", thread_number, state_i);
904 for (
unsigned short i = 0; i < source_functions_size; i++)
906 if (current_state->e[i])
907 if (current_state->e[i]->
used)
909 source_functions_cloned[i]->set_active_element(current_state->e[i]);
910 source_functions_cloned[i]->set_transform(current_state->sub_idx[i]);
918 for (
unsigned short edge = 0; edge < current_state->rep->
nvert; edge++)
921 this->info(
"Thread %i, state %i, edge %i", thread_number, state_i, edge);
923 if (!assemble_everywhere)
925 bool target_marker =
false;
926 for (
unsigned short i = 0; i < internal_markers.size(); i++)
928 if (current_state->rep->
en[edge]->
marker == internal_markers[i])
930 target_marker =
true;
938 memset(result_local, 0,
sizeof(Scalar)* this->number_of_integrals);
941 for (
int i = 0; i < this->number_of_integrals; i++)
942 orders[i] = Hermes::Ord(0);
943 for (
int i = 0; i < source_functions_size; i++)
944 func_ord[i] =
new Func<Ord>(source_functions_cloned[i]->get_fn_order());
946 this->order(func_ord, orders);
948 for (
int i = 0; i < this->number_of_integrals; i++)
951 int order_int = order.get_order();
955 int n = init_surface_geometry_points_allocated(refmap, order_int, edge, current_state->rep->
en[edge]->
marker, geometry, jacobian_x_weights);
957 for (
int i = 0; i < source_functions_size; i++)
959 if (current_state->e[i])
961 if (current_state->e[i]->
used)
962 func[i] =
init_fn(source_functions_cloned[i], order_int);
964 func[i] = init_zero_fn<Scalar>(current_state->rep->get_mode(), order_int);
967 func[i] = init_zero_fn<Scalar>(current_state->rep->get_mode(), order_int);
970 this->integral(n, jacobian_x_weights, func, &geometry, result_local);
972 for (
int i = 0; i < source_functions_size; i++)
978 for (
int i = 0; i < this->number_of_integrals; i++)
979 result_thread_local[i] += .5 * result_local[i];
983 this->add_results(result_thread_local, result);
984 free_with_check(result_thread_local);
986 for (
int i = 0; i < source_functions_size; i++)
987 delete source_functions_cloned[i];
988 free_with_check(source_functions_cloned);
989 free_with_check(orders);
990 free_with_check(func_ord);
991 free_with_check(func);
992 free_with_check(result_local);
996 for (
int i = 0; i < num_states; i++)
998 free_with_check(states);
Scalar * calculate(std::vector< std::string > markers)
#define CENTROID_QUAD_X
Centroid of the reference quadrilateral.
Surface integral calculator.
Represents a function defined on a mesh.
Volumetric integral calculator.
Used to pass the instances of Space around.
bool used
array item usage flag
HERMES_API Func< double > * init_fn(PrecalcShapeset *fu, RefMap *rm, const int order)
Init the shape function for the evaluation of the volumetric/surface integral (transformation of valu...
unsigned char nvert
number of vertices (3 or 4)
Scalar * calculate(std::vector< std::string > markers)
int get_num_dofs() const
Returns the number of basis functions contained in the space.
int get_inv_ref_order() const
Returns the increase in the integration order due to the reference map.
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.
Node * en[H2D_MAX_NUMBER_EDGES]
edge node pointers
#define H2D_GET_H_ORDER(encoded_order)
Macros for combining quad horizontal and vertical encoded_orders.
static void untransform(Element *e, double x, double y, double &xi1, double &xi2)
Struct for return type of get_internal_marker().
void set_quad_2d(Quad2D *quad_2d)
::xsd::cxx::tree::string< char, simple_type > string
C++ type corresponding to the string XML Schema built-in type.
Represents a finite element space over a domain.
Represents the reference mapping.
State ** get_states(std::vector< MeshSharedPtr > meshes, unsigned int &states_count)
#define CENTROID_TRI_X
Centroid of the reference triangle.
IntegralCalculator(MeshFunctionSharedPtr< Scalar > source_function, int number_of_integrals)
virtual Scalar * calculate(std::vector< std::string > markers)=0
Represents the solution of a PDE.
HERMES_API unsigned char init_geometry_points_allocated(RefMap **reference_mapping, unsigned short reference_mapping_count, int order, GeomVol< double > &geometry, double *jacobian_x_weights)
virtual void set_active_element(Element *e)