From 2621c658bfd7106f4abf06cb6a1ecd372c40eed5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 31 Dec 2018 09:28:30 -0700 Subject: [PATCH] Clean up some functions in FEFieldFunction. --- .../numerics/fe_field_function.templates.h | 76 ++++++++++++++----- 1 file changed, 56 insertions(+), 20 deletions(-) diff --git a/include/deal.II/numerics/fe_field_function.templates.h b/include/deal.II/numerics/fe_field_function.templates.h index b80d1bffd9..0f3b8f614e 100644 --- a/include/deal.II/numerics/fe_field_function.templates.h +++ b/include/deal.II/numerics/fe_field_function.templates.h @@ -273,10 +273,9 @@ namespace Functions const unsigned int n_cells = compute_point_locations(points, cells, qpoints, maps); - hp::MappingCollection mapping_collection(mapping); - const hp::FECollection &fe_collection = dh->get_fe_collection(); - hp::QCollection quadrature_collection; + // Create quadrature collection + hp::QCollection quadrature_collection; for (unsigned int i = 0; i < n_cells; ++i) { // Number of quadrature points on this cell @@ -286,21 +285,25 @@ namespace Functions quadrature_collection.push_back(Quadrature(qpoints[i], ww)); } - // Get a function value object - hp::FEValues fe_v(mapping_collection, - fe_collection, + + // Now gather all the information we need + hp::MappingCollection mapping_collection(mapping); + hp::FEValues fe_v(mapping_collection, + dh->get_fe_collection(), quadrature_collection, update_values); - // Now gather all the information we need for (unsigned int i = 0; i < n_cells; ++i) { AssertThrow(!cells[i]->is_artificial(), VectorTools::ExcPointNotAvailableHere()); + fe_v.reinit(cells[i], i, 0); const unsigned int nq = qpoints[i].size(); + std::vector> vvalues( nq, Vector(this->n_components)); fe_v.get_present_fe_values().get_function_values(data_vector, vvalues); + for (unsigned int q = 0; q < nq; ++q) values[maps[i][q]] = vvalues[q]; } @@ -317,10 +320,17 @@ namespace Functions { Assert(points.size() == values.size(), ExcDimensionMismatch(points.size(), values.size())); + + // Simply forward everything to the vector_value_list() + // function. This requires a temporary object, but everything we + // do here is so expensive that that really doesn't make any + // difference any more. std::vector> vvalues( points.size(), Vector(this->n_components)); + vector_value_list(points, vvalues); + for (unsigned int q = 0; q < points.size(); ++q) values[q] = vvalues[q](component); } @@ -343,10 +353,9 @@ namespace Functions const unsigned int n_cells = compute_point_locations(points, cells, qpoints, maps); - hp::MappingCollection mapping_collection(mapping); - const hp::FECollection &fe_collection = dh->get_fe_collection(); - hp::QCollection quadrature_collection; + // Create quadrature collection + hp::QCollection quadrature_collection; for (unsigned int i = 0; i < n_cells; ++i) { // Number of quadrature points on this cell @@ -356,17 +365,20 @@ namespace Functions quadrature_collection.push_back(Quadrature(qpoints[i], ww)); } - // Get a function value object - hp::FEValues fe_v(mapping_collection, - fe_collection, + + // Now gather all the information we need + hp::MappingCollection mapping_collection(mapping); + hp::FEValues fe_v(mapping_collection, + dh->get_fe_collection(), quadrature_collection, update_gradients); - // Now gather all the information we need for (unsigned int i = 0; i < n_cells; ++i) { AssertThrow(!cells[i]->is_artificial(), VectorTools::ExcPointNotAvailableHere()); + fe_v.reinit(cells[i], i, 0); + const unsigned int nq = qpoints[i].size(); std::vector< std::vector>> @@ -375,6 +387,7 @@ namespace Functions this->n_components)); fe_v.get_present_fe_values().get_function_gradients(data_vector, vgrads); + for (unsigned int q = 0; q < nq; ++q) { const unsigned int s = vgrads[q].size(); @@ -385,6 +398,8 @@ namespace Functions } } + + template void FEFieldFunction::gradient_list( @@ -394,16 +409,24 @@ namespace Functions { Assert(points.size() == values.size(), ExcDimensionMismatch(points.size(), values.size())); + + // Simply forward everything to the vector_gradient_list() + // function. This requires a temporary object, but everything we + // do here is so expensive that that really doesn't make any + // difference any more. std::vector>> vvalues(points.size(), std::vector>( this->n_components)); + vector_gradient_list(points, vvalues); + for (unsigned int q = 0; q < points.size(); ++q) values[q] = vvalues[q][component]; } + template void FEFieldFunction::vector_laplacian_list( @@ -419,10 +442,9 @@ namespace Functions const unsigned int n_cells = compute_point_locations(points, cells, qpoints, maps); - hp::MappingCollection mapping_collection(mapping); - const hp::FECollection &fe_collection = dh->get_fe_collection(); - hp::QCollection quadrature_collection; + // Create quadrature collection + hp::QCollection quadrature_collection; for (unsigned int i = 0; i < n_cells; ++i) { // Number of quadrature points on this cell @@ -432,9 +454,11 @@ namespace Functions quadrature_collection.push_back(Quadrature(qpoints[i], ww)); } - // Get a function value object - hp::FEValues fe_v(mapping_collection, - fe_collection, + + // Now gather all the information we need + hp::MappingCollection mapping_collection(mapping); + hp::FEValues fe_v(mapping_collection, + dh->get_fe_collection(), quadrature_collection, update_hessians); // Now gather all the information we need @@ -442,17 +466,22 @@ namespace Functions { AssertThrow(!cells[i]->is_artificial(), VectorTools::ExcPointNotAvailableHere()); + fe_v.reinit(cells[i], i, 0); + const unsigned int nq = qpoints[i].size(); std::vector> vvalues( nq, Vector(this->n_components)); fe_v.get_present_fe_values().get_function_laplacians(data_vector, vvalues); + for (unsigned int q = 0; q < nq; ++q) values[maps[i][q]] = vvalues[q]; } } + + template void FEFieldFunction::laplacian_list( @@ -462,10 +491,17 @@ namespace Functions { Assert(points.size() == values.size(), ExcDimensionMismatch(points.size(), values.size())); + + // Simply forward everything to the vector_gradient_list() + // function. This requires a temporary object, but everything we + // do here is so expensive that that really doesn't make any + // difference any more. std::vector> vvalues( points.size(), Vector(this->n_components)); + vector_laplacian_list(points, vvalues); + for (unsigned int q = 0; q < points.size(); ++q) values[q] = vvalues[q](component); } -- 2.39.5