From ad1ee5e7fbabc1abe7b0bb8828e830b526212bc4 Mon Sep 17 00:00:00 2001 From: Mikael Mortensen Date: Fri, 11 Oct 2019 10:25:03 +0200 Subject: [PATCH] Compatible with fenics 2019.1 --- fenicstools/fem/cr_divergence.cpp | 56 ++++++++++--------- fenicstools/fem/gradient_weight.cpp | 83 ++++++++++++++--------------- fenicstools/fem/interpolation.cpp | 66 +++++++++++------------ tests/test_Probe.py | 5 ++ 4 files changed, 105 insertions(+), 105 deletions(-) diff --git a/fenicstools/fem/cr_divergence.cpp b/fenicstools/fem/cr_divergence.cpp index ad051c7..60050f1 100644 --- a/fenicstools/fem/cr_divergence.cpp +++ b/fenicstools/fem/cr_divergence.cpp @@ -24,7 +24,7 @@ cfg['library_dirs'] = dolfin_pc['library_dirs'] using namespace dolfin; -// Base case for all divergence computations. +// Base case for all divergence computations. // Compute divergence of vector field u. void cr_divergence(Function& divu, const Function& u) { @@ -36,7 +36,7 @@ void cr_divergence(Function& divu, const Function& u) std::shared_ptr U = u.vector(); std::shared_ptr DIVU = divu.vector(); - // Figure out about the local dofs of DG0 + // Figure out about the local dofs of DG0 std::pair first_last_dof = DG0_dofmap->ownership_range(); std::size_t first_dof = first_last_dof.first; @@ -48,7 +48,7 @@ void cr_divergence(Function& divu, const Function& u) // Get topological dimension so that we know what Facet is const Mesh mesh = *u.function_space()->mesh(); - std::size_t tdim = mesh.topology().dim(); + std::size_t tdim = mesh.topology().dim(); // Get the info on length of u and gdim for the dot product std::size_t gdim = mesh.geometry().dim(); std::size_t udim = u.value_dimension(0); // Rank 1! @@ -65,10 +65,10 @@ void cr_divergence(Function& divu, const Function& u) { Point cell_mp = cell->midpoint(); double cell_volume = cell->volume(); - + // Dofs of CR on all facets of the cell, global order auto facets_dofs = CR1_dofmap->cell_dofs(cell->index()); - + double cell_integral = 0; std::size_t local_facet_index = 0; for(FacetIterator facet(*cell); !facet.end(); ++facet) @@ -79,7 +79,7 @@ void cr_divergence(Function& divu, const Function& u) else if(tdim == 3) facet_measure = Face(mesh, facet->index()).area(); // Tdim 1 will not happen because CR is not defined there - + Point facet_normal = facet->normal(); // Flip the normal if it is not outer already @@ -115,7 +115,7 @@ void cr_divergence(Function& divu, const Function& u) void cr_divergence2(Function& divu, const Function& u, const FunctionSpace& DGscalar, const FunctionSpace& CRvector) { - // Divu scalar from u vector + // Divu scalar from u vector std::size_t u_rank = u.value_rank(); std::size_t divu_rank = divu.value_rank(); if((divu_rank == 0) and (u_rank == 1)) @@ -139,10 +139,10 @@ void cr_divergence2(Function& divu, const Function& u, std::shared_ptr U = u.vector(); std::shared_ptr DIVU = divu.vector(); - // Get dofmaps + // Get dofmaps std::shared_ptr DGvector = divu.function_space(); std::shared_ptr CRtensor = u.function_space(); - + std::size_t len_divu = divu.value_dimension(0); // With scalar U can be extracted only once std::vector U_values; @@ -163,7 +163,7 @@ void cr_divergence2(Function& divu, const Function& u, CRi_rows = CRvector[i]->dofmap()->dofs(); std::size_t m = CRi_rows.size(); U_i->zero(); - U_i->set(U_values.data(), m, CRi_rows.data()); + U_i->set(U_values.data(), m, CRi_rows.data()); U_i->apply("insert"); } // For tensor, U_i represents T_ij, (T_i0, T_i1, T_i2) @@ -183,7 +183,7 @@ void cr_divergence2(Function& divu, const Function& u, std::vector CRj_rows = CRvector[j]->dofmap()->dofs(); - U_i->set(U_values.data(), m, CRj_rows.data()); + U_i->set(U_values.data(), m, CRj_rows.data()); U_i->apply("insert"); } } @@ -197,22 +197,22 @@ void cr_divergence2(Function& divu, const Function& u, std::size_t m = DGi_rows.size(); std::vector DIVU_i_values; DIVU_i->get_local(DIVU_i_values); - DIVU->set(DIVU_i_values.data(), m, DGi_rows.data()); + DIVU->set(DIVU_i_values.data(), m, DGi_rows.data()); DIVU->apply("insert"); } } } -// Base case for all divergence computations. +// Base case for all divergence computations. // Compute divergence of vector field u. void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A, - const FunctionSpace& DGscalar, + const FunctionSpace& DGscalar, const FunctionSpace& CRvector) { std::shared_ptr CR1_dofmap = CRvector.dofmap(), DG0_dofmap = DGscalar.dofmap(); - // Figure out about the local dofs of DG0 + // Figure out about the local dofs of DG0 std::pair first_last_dof = DG0_dofmap->ownership_range(); std::size_t first_dof = first_last_dof.first; @@ -221,25 +221,25 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A, // Get topological dimension so that we know what Facet is const Mesh mesh = *DGscalar.mesh(); - std::size_t tdim = mesh.topology().dim(); + std::size_t tdim = mesh.topology().dim(); std::size_t gdim = mesh.geometry().dim(); - + std::vector columns; std::vector values; - + // Fill the values for(CellIterator cell(mesh); !cell.end(); ++cell) { auto dg_dofs = DG0_dofmap->cell_dofs(cell->index()); // There is only one DG0 dof per cell dolfin::la_index cell_dof = dg_dofs[0]; - + Point cell_mp = cell->midpoint(); double cell_volume = cell->volume(); - std::size_t local_facet_index = 0; - + std::size_t local_facet_index = 0; + auto cr_dofs = CR1_dofmap->cell_dofs(cell->index()); - + for(FacetIterator facet(*cell); !facet.end(); ++facet) { double facet_measure=0; @@ -248,23 +248,23 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A, else if(tdim == 3) facet_measure = Face(mesh, facet->index()).area(); // Tdim 1 will not happen because CR is not defined there - + Point facet_normal = facet->normal(); // Flip the normal if it is not outer already Point facet_mp = facet->midpoint(); double sign = (facet_normal.dot(facet_mp - cell_mp) > 0.0) ? 1.0 : -1.0; facet_normal *= (sign*facet_measure/cell_volume); - + // Dofs of CR on the facet, local order std::vector facet_dofs; CR1_dofmap->tabulate_facet_dofs(facet_dofs, local_facet_index); - + for (std::size_t j = 0 ; j < facet_dofs.size(); j++) - { + { columns.push_back(cr_dofs[facet_dofs[j]]); values.push_back(facet_normal[j]); - } + } local_facet_index += 1; } M.setrow(cell_dof, columns, values); @@ -275,7 +275,6 @@ void cr_divergence_matrix(GenericMatrix& M, GenericMatrix& A, //std::shared_ptr Cp = MatMatMult(M, A); //return Cp; } - namespace py = pybind11; PYBIND11_MODULE(cr_divergence, m) @@ -299,4 +298,3 @@ PYBIND11_MODULE(cr_divergence, m) cr_divergence_matrix(M, A, _U, _V); }); } - diff --git a/fenicstools/fem/gradient_weight.cpp b/fenicstools/fem/gradient_weight.cpp index e38d744..37762f4 100644 --- a/fenicstools/fem/gradient_weight.cpp +++ b/fenicstools/fem/gradient_weight.cpp @@ -25,27 +25,27 @@ using namespace dolfin; void compute_weight(Function& DG) { // Compute weights for averaging with neighboring cells - + // Get the mesh, element and dofmap - std::shared_ptr V = DG.function_space(); + std::shared_ptr V = DG.function_space(); std::shared_ptr mesh = V->mesh(); std::shared_ptr element = V->element(); std::shared_ptr dofmap_u = V->dofmap(); - + // Allocate storage for weights on one cell - std::vector ws(element->space_dimension()); - + std::vector ws(element->space_dimension()); + // Compute weights - GenericVector& dg_vector = *DG.vector(); + GenericVector& dg_vector = *DG.vector(); dg_vector.zero(); for (CellIterator cell(*mesh, "all"); !cell.end(); ++cell) { auto dofs = dofmap_u->cell_dofs(cell->index()); - + std::fill(ws.begin(), ws.end(), 1./cell->volume()); - dg_vector.add_local(ws.data(), dofs.size(), dofs.data()); - } - dg_vector.apply("insert"); + dg_vector.add_local(ws.data(), dofs.size(), dofs.data()); + } + dg_vector.apply("insert"); } std::size_t dof_owner(std::vector > all_ranges, @@ -66,7 +66,7 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG) std::vector values; std::vector > allcolumns; std::vector > allvalues; - + const std::pair row_range = A.local_range(0); const std::size_t m = row_range.second - row_range.first; GenericVector& weight = *DG.vector(); @@ -75,26 +75,26 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG) weight_range_vec[0] = weight_range.first; weight_range_vec[1] = weight_range.second; int dm = weight_range.second-weight_range.first; - + const MPI_Comm mpi_comm = DG.function_space()->mesh()->mpi_comm(); - + // Communicate local_ranges of weights std::vector > all_ranges; dolfin::MPI::all_gather(mpi_comm, weight_range_vec, all_ranges); - + // Number of MPI processes std::size_t num_processes = dolfin::MPI::size(mpi_comm); // Some weights live on other processes and need to be communicated // Create list of off-process weights - std::vector > dofs_needed(num_processes); + std::vector > dofs_needed(num_processes); for (std::size_t row = 0; row < m; row++) - { + { // Get global row number const std::size_t global_row = row + row_range.first; - + A.getrow(global_row, columns, values); - + for (std::size_t i = 0; i < columns.size(); i++) { std::size_t dof = columns[i]; @@ -109,14 +109,14 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG) // Communicate to all which weights are needed by the process std::vector > dofs_needed_recv; dolfin::MPI::all_to_all(mpi_comm, dofs_needed, dofs_needed_recv); - + // Fetch the weights that must be communicated - std::vector > weights_to_send(num_processes); + std::vector > weights_to_send(num_processes); for (std::size_t p = 0; p < num_processes; p++) { - if (p == dolfin::MPI::rank(mpi_comm)) + if (p == dolfin::MPI::rank(mpi_comm)) continue; - + std::vector dofs = dofs_needed_recv[p]; std::map send_weights; for (std::size_t k = 0; k < dofs.size(); k++) @@ -126,25 +126,25 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG) } std::vector > weights_to_send_recv; dolfin::MPI::all_to_all(mpi_comm, weights_to_send, weights_to_send_recv); - + // Create a map for looking up received weights std::map received_weights; for (std::size_t p = 0; p < num_processes; p++) { if (p == dolfin::MPI::rank(mpi_comm)) continue; - + for (std::size_t k = 0; k < dofs_needed[p].size(); k++) { - received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k]; + received_weights[dofs_needed[p][k]] = weights_to_send_recv[p][k]; } } - + for (std::size_t row = 0; row < m; row++) - { + { // Get global row number const std::size_t global_row = row + row_range.first; - + A.getrow(global_row, columns, values); for (std::size_t i = 0; i < values.size(); i++) { @@ -155,14 +155,14 @@ void compute_DG0_to_CG_weight_matrix(GenericMatrix& A, Function& DG) } else { - values[i] = weight[columns[i]-weight_range.first]; + values[i] = weight[columns[i]-weight_range.first]; } // values[i] = 1./values[i]; } - + double s = std::accumulate(values.begin(), values.end(), 0.0); std::transform(values.begin(), values.end(), values.begin(), - std::bind2nd(std::multiplies(), 1./s)); + std::bind2nd(std::multiplies(), 1./s)); for (std::size_t i=0; i MatTransposeMatMult(const GenericMatrix& A, const GenericMatrix& B) { @@ -210,7 +210,7 @@ std::shared_ptr MatMatTransposeMult(const GenericMatrix& A, const const PETScMatrix* Bp = &as_type(B); Mat CC; PetscErrorCode ierr = MatMatTransposeMult(Ap->mat(), Bp->mat(), MAT_INITIAL_MATRIX, PETSC_DEFAULT, &CC); - return PETScMatrix(CC).copy(); + return PETScMatrix(CC).copy(); } // std::shared_ptr matMatMult(const GenericMatrix& A, const GenericMatrix& B) @@ -248,4 +248,3 @@ PYBIND11_MODULE(gradient_weight, m) return MatTransposeMatMult(A, B); }); } - diff --git a/fenicstools/fem/interpolation.cpp b/fenicstools/fem/interpolation.cpp index 42b64ce..78709e1 100644 --- a/fenicstools/fem/interpolation.cpp +++ b/fenicstools/fem/interpolation.cpp @@ -11,7 +11,6 @@ cfg['compiler_args'] = ['-std=c++11', '-DHAS_MPI'] #include #include #include - #include #include #include @@ -54,14 +53,14 @@ struct lt_coordinate } // Tolerance - const double TOL; + const double TOL; }; -void extract_dof_component_map(std::unordered_map& dof_component_map, - const FunctionSpace& V, +void extract_dof_component_map(std::unordered_map& dof_component_map, + const FunctionSpace& V, int* component) -{ +{ // Extract sub dofmaps recursively and store dof to component map if (V.element()->num_sub_elements() == 0) { @@ -101,10 +100,10 @@ bool in_bounding_box(const std::vector& point, } } return true;} - -std::map, std::vector, lt_coordinate> + +std::map, std::vector, lt_coordinate> tabulate_coordinates_to_dofs(const FunctionSpace& V) -{ +{ std::map, std::vector, lt_coordinate> coords_to_dofs(lt_coordinate(1.0e-12)); @@ -166,8 +165,8 @@ tabulate_coordinates_to_dofs(const FunctionSpace& V) return coords_to_dofs; } -void interpolate1(const Expression& u0, Function& u) -{ +void interpolate1(const Expression& u0, Function& u) +{ // Get function space and element interpolating to dolfin_assert(u.function_space()); const FunctionSpace& V = *u.function_space(); @@ -237,7 +236,7 @@ void interpolate1(const Expression& u0, Function& u) u.vector()->apply("insert"); } -void interpolate2(const Function& u0, Function& u) +void interpolate2(const Function& u0, Function& u) { // Interpolate from Function u0 to Function u. // This mesh of u0 may be different from that of u @@ -473,7 +472,7 @@ void interpolate_any1(const Function& u0, Function& u) // be required // 6) Call evaluate_dofs for all cells on local mesh using the // wrapped function as the ufc function. - // 7) Update coefficients of local vector with results from + // 7) Update coefficients of local vector with results from // evaluate_dofs // Get function spaces of Functions interpolating to/from @@ -550,7 +549,7 @@ void interpolate_any1(const Function& u0, Function& u) // Create map from coordinate to global result of eval static std::map, std::vector, lt_coordinate> coords_to_values(lt_coordinate(1.0e-12)); - + // Create an Expression wrapping the set of coordinates // Objective is to obtain all points in need of evaluation. static std::set> coords; @@ -571,7 +570,7 @@ void interpolate_any1(const Function& u0, Function& u) cord._xx = x; ufc::cell ufc_cell; std::vector vertex_coordinates; - std::vector cell_coefficients(dofmap.max_element_dofs()); + std::vector cell_coefficients(dofmap.max_element_dofs()); const std::size_t local_size = dofmap.ownership_range().second - dofmap.ownership_range().first; for (CellIterator cell(mesh1); !cell.end(); ++cell) @@ -589,8 +588,8 @@ void interpolate_any1(const Function& u0, Function& u) // for (auto y = it->begin(); y != it->end(); ++y) // std::cout << *y << " "; // std::cout << std::endl; - // } - // std::cout << coords.size() << " " << coords_to_dofs.size() << std::endl; + // } + // std::cout << coords.size() << " " << coords_to_dofs.size() << std::endl; // Search this process first for all coordinates in u's local mesh std::vector points_not_found; @@ -698,22 +697,22 @@ void interpolate_any1(const Function& u0, Function& u) coords_to_values.insert(std::make_pair(x, values)); } } - + // Create an Expression wrapping coords_to_values such that these // results will be retrieved when calling evaluate_dofs class WrapperFunction : public Expression { public: - + WrapperFunction(int value_shape) : Expression(value_shape) {}; - + mutable std::vector _xx; void eval(Array& values, const Array& x) const { for (uint j = 0; j < x.size(); j++) _xx[j] = x[j]; - + const std::vector& v = coords_to_values[_xx]; for (std::size_t j = 0; j < v.size(); j++) values[j] = v[j]; @@ -722,12 +721,12 @@ void interpolate_any1(const Function& u0, Function& u) WrapperFunction wrapper(u.value_size()); wrapper._xx = x; - + // Iterate over mesh and interpolate on each cell // ufc::cell ufc_cell; // std::vector vertex_coordinates; // std::vector cell_coefficients(dofmap.max_element_dofs()); - // + // // const std::size_t local_size = dofmap.ownership_range().second // - dofmap.ownership_range().first; for (CellIterator cell(mesh1); !cell.end(); ++cell) @@ -736,7 +735,7 @@ void interpolate_any1(const Function& u0, Function& u) cell->get_vertex_coordinates(vertex_coordinates); cell->get_cell_data(ufc_cell); - // Call evaluate_dofs with wrapper function around the globally + // Call evaluate_dofs with wrapper function around the globally // computed interpolation points. element.evaluate_dofs(cell_coefficients.data(), wrapper, vertex_coordinates.data(), ufc_cell.orientation, ufc_cell); @@ -752,11 +751,11 @@ void interpolate_any1(const Function& u0, Function& u) local_u_vector[d] = cell_coefficients[i]; } } - + // Set and finalize vector u.vector()->set_local(local_u_vector); u.vector()->apply("insert"); -} +} // void interpolate_any2(const Expression& u0, Function& u) { @@ -811,10 +810,10 @@ void interpolate_any2(const Expression& u0, Function& u) // Get dofmap of u dolfin_assert(V1.dofmap()); const GenericDofMap& dofmap = *V1.dofmap(); - + ufc::cell ufc_cell; std::vector vertex_coordinates; - std::vector cell_coefficients(dofmap.max_element_dofs()); + std::vector cell_coefficients(dofmap.max_element_dofs()); const std::size_t local_size = dofmap.ownership_range().second - dofmap.ownership_range().first; @@ -826,7 +825,7 @@ void interpolate_any2(const Expression& u0, Function& u) cell->get_vertex_coordinates(vertex_coordinates); cell->get_cell_data(ufc_cell); - // Call evaluate_dofs with wrapper function around the globally + // Call evaluate_dofs with wrapper function around the globally // computed interpolation points. element.evaluate_dofs(cell_coefficients.data(), u0, vertex_coordinates.data(), ufc_cell.orientation, ufc_cell); @@ -842,17 +841,17 @@ void interpolate_any2(const Expression& u0, Function& u) local_u_vector[d] = cell_coefficients[i]; } } - + // Set and finalize vector u.vector()->set_local(local_u_vector); u.vector()->apply("insert"); -} +} PYBIND11_MODULE(interpolation, m) { - m.def("interpolate", (void (*)(const Function&, Function&)) + m.def("interpolate", (void (*)(const Function&, Function&)) &interpolate2); - m.def("interpolate", (void (*)(const Expression&, Function&)) + m.def("interpolate", (void (*)(const Expression&, Function&)) &interpolate1); m.def("interpolate", [](py::object u0, py::object v){ auto _u = u0.attr("_cpp_object"); @@ -883,4 +882,3 @@ PYBIND11_MODULE(interpolation, m) } }); } - diff --git a/tests/test_Probe.py b/tests/test_Probe.py index 3af16ef..0116d2c 100755 --- a/tests/test_Probe.py +++ b/tests/test_Probe.py @@ -51,3 +51,8 @@ def test_Probe_vectorfunctionspace_3D(VF3_self): assert round(p[0][0] - 0.25, 7) == 0 assert round(p[1][1] - 0.50, 7) == 0 assert round(p[1][2] - 0.75, 7) == 0 + +if __name__ == '__main__': + mesh = UnitSquareMesh(4, 4) + V = FunctionSpace(mesh, 'CG', 1) + test_Probe_functionspace_2D(V) \ No newline at end of file