Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 : #include "libmesh/libmesh_common.h"
19 :
20 : #ifdef LIBMESH_HAVE_PETSC
21 :
22 : // Local Includes
23 : #include "libmesh/petsc_preconditioner.h"
24 : #include "libmesh/petsc_macro.h"
25 : #include "libmesh/petsc_matrix.h"
26 : #include "libmesh/petsc_vector.h"
27 : #include "libmesh/libmesh_common.h"
28 : #include "libmesh/enum_preconditioner_type.h"
29 : #include "libmesh/elem.h"
30 : #include "libmesh/equation_systems.h"
31 : #include "libmesh/dof_map.h"
32 :
33 : namespace libMesh
34 : {
35 :
36 : template <typename T>
37 0 : PetscPreconditioner<T>::PetscPreconditioner (const libMesh::Parallel::Communicator & comm_in) :
38 0 : Preconditioner<T>(comm_in)
39 0 : {}
40 :
41 :
42 :
43 : template <typename T>
44 0 : void PetscPreconditioner<T>::apply(const NumericVector<T> & x, NumericVector<T> & y)
45 : {
46 0 : PetscVector<T> & x_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(x));
47 0 : PetscVector<T> & y_pvec = cast_ref<PetscVector<T> &>(const_cast<NumericVector<T> &>(y));
48 :
49 0 : Vec x_vec = x_pvec.vec();
50 0 : Vec y_vec = y_pvec.vec();
51 :
52 0 : LibmeshPetscCall(PCApply(_pc, x_vec, y_vec));
53 0 : }
54 :
55 :
56 :
57 :
58 : template <typename T>
59 0 : void PetscPreconditioner<T>::init ()
60 : {
61 0 : libmesh_error_msg_if(!this->_matrix, "ERROR: No matrix set for PetscPreconditioner, but init() called");
62 :
63 : // Clear the preconditioner in case it has been created in the past
64 0 : if (!this->_is_initialized)
65 : {
66 : // Should probably use PCReset(), but it's not working at the moment so we'll destroy instead
67 0 : if (_pc)
68 0 : _pc.destroy();
69 :
70 0 : LibmeshPetscCall(PCCreate(this->comm().get(), _pc.get()));
71 :
72 0 : auto pmatrix = cast_ptr<PetscMatrixBase<T> *>(this->_matrix);
73 0 : _mat = pmatrix->mat();
74 : }
75 :
76 0 : LibmeshPetscCall(PCSetOperators(_pc, _mat, _mat));
77 :
78 : // Set the PCType. Note: this used to be done *before* the call to
79 : // PCSetOperators(), and only when !_is_initialized, but
80 : // 1.) Some preconditioners (those employing sub-preconditioners,
81 : // for example) have to call PCSetUp(), and can only do this after
82 : // the operators have been set.
83 : // 2.) It should be safe to call set_petsc_preconditioner_type()
84 : // multiple times.
85 0 : set_petsc_preconditioner_type(this->_preconditioner_type, *_pc);
86 :
87 0 : this->_is_initialized = true;
88 0 : }
89 :
90 :
91 :
92 : template <typename T>
93 0 : void PetscPreconditioner<T>::clear()
94 : {
95 : // Calls custom deleter
96 0 : _pc.destroy();
97 0 : }
98 :
99 :
100 :
101 : template <typename T>
102 0 : PC PetscPreconditioner<T>::pc()
103 : {
104 0 : return _pc;
105 : }
106 :
107 :
108 :
109 : template <typename T>
110 55429 : void PetscPreconditioner<T>::set_petsc_preconditioner_type (const PreconditionerType & preconditioner_type, PC & pc)
111 : {
112 : // get the communicator from the PETSc object
113 : Parallel::communicator comm;
114 55429 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, & comm);
115 55429 : if (ierr != LIBMESH_PETSC_SUCCESS)
116 0 : libmesh_error_msg("Error retrieving communicator");
117 :
118 : #define CasePCSetType(PreconditionerType, PCType) \
119 : case PreconditionerType: \
120 : LibmeshPetscCallA(comm, PCSetType (pc, const_cast<KSPType>(PCType))); \
121 : break;
122 :
123 55429 : switch (preconditioner_type)
124 : {
125 140 : CasePCSetType(IDENTITY_PRECOND, PCNONE)
126 0 : CasePCSetType(CHOLESKY_PRECOND, PCCHOLESKY)
127 0 : CasePCSetType(ICC_PRECOND, PCICC)
128 733 : CasePCSetType(ILU_PRECOND, PCILU)
129 0 : CasePCSetType(LU_PRECOND, PCLU)
130 0 : CasePCSetType(ASM_PRECOND, PCASM)
131 0 : CasePCSetType(JACOBI_PRECOND, PCJACOBI)
132 50706 : CasePCSetType(BLOCK_JACOBI_PRECOND, PCBJACOBI)
133 0 : CasePCSetType(SOR_PRECOND, PCSOR)
134 0 : CasePCSetType(EISENSTAT_PRECOND, PCEISENSTAT)
135 0 : CasePCSetType(AMG_PRECOND, PCHYPRE)
136 0 : CasePCSetType(SVD_PRECOND, PCSVD)
137 0 : CasePCSetType(USER_PRECOND, PCMAT)
138 3850 : CasePCSetType(SHELL_PRECOND, PCSHELL)
139 :
140 0 : default:
141 0 : libMesh::err << "ERROR: Unsupported PETSC Preconditioner: "
142 0 : << Utility::enum_to_string(preconditioner_type) << std::endl
143 0 : << "Continuing with PETSC defaults" << std::endl;
144 : }
145 :
146 : // Set additional options if we are doing AMG and
147 : // HYPRE is available
148 : #ifdef LIBMESH_HAVE_PETSC_HYPRE
149 55429 : if (preconditioner_type == AMG_PRECOND)
150 0 : LibmeshPetscCallA(comm, PCHYPRESetType(pc, "boomeramg"));
151 : #endif
152 :
153 : // Let the commandline override stuff
154 55429 : LibmeshPetscCallA(comm, PCSetFromOptions(pc));
155 55429 : }
156 :
157 :
158 :
159 : #ifdef LIBMESH_HAVE_PETSC_HYPRE
160 : template <typename T>
161 259204 : void PetscPreconditioner<T>::set_petsc_aux_data(PC & pc, System & sys, const unsigned v)
162 : {
163 : // Get the communicator from the PETSc object
164 : Parallel::communicator comm;
165 259204 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
166 259204 : libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
167 : "Error retrieving communicator");
168 :
169 : // Make sure the preconditioner options are set
170 259204 : LibmeshPetscCallA(comm, PCSetFromOptions(pc));
171 :
172 : // Get the type of preconditioner we are using
173 259204 : PCType pc_type = nullptr;
174 259204 : LibmeshPetscCallA(comm, PCGetType(pc, &pc_type));
175 :
176 : // Check if hypre ams/ads, otherwise we quit with nothing to do
177 1012936 : if (pc_type && std::string(pc_type) == PCHYPRE)
178 : {
179 : // Get the hypre preconditioner we are using
180 804 : PCType hypre_type = nullptr;
181 804 : LibmeshPetscCallA(comm, PCHYPREGetType(pc, &hypre_type));
182 :
183 : // If not ams/ads, we quit with nothing to do
184 1608 : if (std::string(hypre_type) == "ams")
185 : {
186 : // If multiple variables, we error out as senseless
187 800 : libmesh_error_msg_if(sys.n_vars() > 1,
188 : "Error applying hypre AMS to a system with multiple "
189 : "variables");
190 : // If not a 1st order Nédélec or a 2d 1st order Raviart-Thomas system, we
191 : // error out as we do not support anything else at the moment
192 816 : libmesh_error_msg_if(sys.variable(v).type() != FEType(1, NEDELEC_ONE) &&
193 : (sys.variable(v).type() != FEType(1, RAVIART_THOMAS) ||
194 : sys.get_mesh().mesh_dimension() != 2),
195 : "Error applying hypre AMS to a system "
196 : "whose variable is not 1st order Nedelec or 1st "
197 : "order Raviart-Thomas on a 2d mesh");
198 800 : set_hypre_ams_data(pc, sys, v);
199 : }
200 8 : else if (std::string(hypre_type) == "ads")
201 : {
202 : // If multiple variables, we error out as senseless
203 4 : libmesh_error_msg_if(sys.n_vars() > 1,
204 : "Error applying hypre ADS to a system with multiple "
205 : "variables");
206 : // If not a 3d 1st order Raviart-Thomas system, we error out as we do not
207 : // support anything else at the moment
208 8 : libmesh_error_msg_if(sys.variable(v).type() != FEType(1, RAVIART_THOMAS) ||
209 : sys.get_mesh().mesh_dimension() != 3,
210 : "Error applying hypre ADS to a system "
211 : "whose variable is not 1st "
212 : "order Raviart-Thomas on a 3d mesh");
213 4 : set_hypre_ads_data(pc, sys, v);
214 : }
215 : }
216 259204 : }
217 :
218 :
219 :
220 : template <typename T>
221 800 : void PetscPreconditioner<T>::set_hypre_ams_data(PC & pc, System & sys, const unsigned v)
222 : {
223 : // Get the communicator from the PETSc object
224 : Parallel::communicator comm;
225 800 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
226 800 : libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
227 : "Error retrieving communicator");
228 0 : Parallel::Communicator Comm(comm);
229 :
230 : // The mesh/problem dimension
231 800 : const unsigned dim = sys.get_mesh().mesh_dimension();
232 :
233 : // Dummy Lagrange system defined over the same mesh so we can enumerate the vertices
234 800 : System & lagrange_sys = sys.get_equation_systems().add_system<System>("__hypre_ams_vertices");
235 800 : lagrange_sys.hide_output() = true;
236 800 : lagrange_sys.add_variable("__lagrange");
237 800 : lagrange_sys.reinit_mesh();
238 :
239 : // Global (i.e. total) and local (i.e. to this processor) number of edges and vertices
240 800 : const std::vector<dof_id_type> n_all_edges = sys.get_dof_map().n_dofs_per_processor(v);
241 800 : const dof_id_type n_glb_edges = std::accumulate(n_all_edges.begin(), n_all_edges.end(), 0);
242 800 : const dof_id_type n_loc_edges = n_all_edges[global_processor_id()];
243 800 : const dof_id_type n_glb_verts = lagrange_sys.n_dofs();
244 800 : const dof_id_type n_loc_verts = lagrange_sys.n_local_dofs();
245 :
246 : // We require local indexing through the vertices and contiguous indexing
247 : // through the edges: since the vertices are enumerated on their own system,
248 : // a local index can simply be obtained by subtracting those vertices in all
249 : // preceding processors; for the edges, if the system we are given by the
250 : // user has multiple variables/splits, we effectively need to add local edge
251 : // numbers to those edges in all preceding processors
252 :
253 : // The number of vertices and edges in all preceding processors
254 0 : dof_id_signed_type vert_offset = -lagrange_sys.get_dof_map().first_dof();
255 800 : dof_id_signed_type edge_offset =
256 0 : std::accumulate(n_all_edges.begin(), n_all_edges.begin() + global_processor_id(), 0);
257 :
258 : // Whether dofs are in variable-major or node-major order
259 800 : bool var_major = !libMesh::on_command_line("--node-major-dofs");
260 :
261 : // If variable-major order, we need only subtract all the preceding dofs; but
262 : // if node-major order, we need to enumerate the dofs on this processor
263 0 : std::vector<dof_id_type> idx_edges;
264 800 : if (var_major)
265 : {
266 800 : edge_offset -= sys.get_dof_map().first_dof();
267 800 : for (auto i : make_range(v))
268 0 : edge_offset -= sys.get_dof_map().n_local_dofs(i);
269 : }
270 : else
271 0 : sys.get_dof_map().local_variable_indices(idx_edges, sys.get_mesh(), v);
272 :
273 : // Create the discrete grandient matrix, representing the edges in terms of its vertices
274 : // Preallocate 2 diagonal + 2 off-diagonal nonzeros as the vertices could fall on either
275 800 : PetscMatrix<Real> G(Comm, n_glb_edges, n_glb_verts, n_loc_edges, n_loc_verts, 2, 2);
276 :
277 : // Create array for the coordinates of the local vertices
278 : PetscReal * coords;
279 800 : LibmeshPetscCallA(comm, PetscMalloc1(dim * n_loc_verts, &coords));
280 :
281 : // Populate the discrete gradient matrix and the coordinates array
282 149100 : for (const auto & elem : sys.get_mesh().active_local_element_ptr_range())
283 401562 : for (auto edge : make_range(elem->n_edges()))
284 : {
285 : // The edge's two vertices
286 : dof_id_type vert_dofs[2];
287 984636 : for (auto vert : make_range(2))
288 : {
289 656424 : const unsigned loc_vert_node_id = elem->local_edge_node(edge, vert);
290 0 : libmesh_assert(elem->is_vertex(loc_vert_node_id));
291 :
292 656424 : const Node & vert_node = elem->node_ref(loc_vert_node_id);
293 656424 : vert_dofs[vert] = vert_node.dof_number(lagrange_sys.number(), 0, 0);
294 :
295 : // If owned, populate coordinates array
296 656424 : if (vert_node.processor_id() == global_processor_id())
297 : {
298 0 : const dof_id_type loc_vert_dof = vert_offset + vert_dofs[vert];
299 0 : libmesh_assert(loc_vert_dof < n_loc_verts);
300 :
301 1902888 : for (auto d : make_range(dim))
302 1360376 : coords[dim * loc_vert_dof + d] = vert_node(d);
303 : }
304 : }
305 :
306 : // The edge's (middle) node
307 328212 : const unsigned loc_edge_node_id = elem->local_edge_node(edge, 2);
308 0 : libmesh_assert(elem->is_edge(loc_edge_node_id));
309 :
310 328212 : const Node & edge_node = elem->node_ref(loc_edge_node_id);
311 :
312 : // If owned, populate discrete gradient matrix
313 328212 : if (edge_node.processor_id() == global_processor_id())
314 : {
315 303502 : const dof_id_type edge_dof = edge_node.dof_number(sys.number(), v, 0);
316 :
317 303502 : const dof_id_type cont_edge_dof = edge_offset +
318 303502 : (var_major ? edge_dof
319 0 : : std::distance(idx_edges.begin(),
320 : std::find(idx_edges.begin(), idx_edges.end(), edge_dof)));
321 :
322 303502 : const Real sign = elem->positive_edge_orientation(edge) ? 1 : -1;
323 :
324 0 : libmesh_assert(cont_edge_dof >= G.row_start());
325 0 : libmesh_assert(cont_edge_dof < G.row_stop());
326 303502 : G.set(cont_edge_dof, vert_dofs[0], sign);
327 303502 : G.set(cont_edge_dof, vert_dofs[1], -sign);
328 : }
329 : }
330 :
331 : // Assemble the discrete gradient matrix
332 800 : G.close();
333 :
334 : // Hand over the matrix and coordinates array
335 800 : LibmeshPetscCallA(comm, PCHYPRESetDiscreteGradient(pc, G.mat()));
336 800 : LibmeshPetscCallA(comm, PCSetCoordinates(pc, dim, n_loc_verts, coords));
337 :
338 : // Free allocated memory for the coordinates array
339 800 : LibmeshPetscCallA(comm, PetscFree(coords));
340 1600 : }
341 :
342 :
343 :
344 : template <typename T>
345 4 : void PetscPreconditioner<T>::set_hypre_ads_data(PC & pc, System & sys, const unsigned v)
346 : {
347 : // Get the communicator from the PETSc object
348 : Parallel::communicator comm;
349 4 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
350 4 : libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
351 : "Error retrieving communicator");
352 0 : Parallel::Communicator Comm(comm);
353 :
354 : // The mesh/problem dimension
355 4 : const unsigned dim = sys.get_mesh().mesh_dimension();
356 :
357 : // Dummy Lagrange and Nédélec systems defined over the same mesh so we can
358 : // enumerate the vertices and edges, respectively
359 4 : System & lagrange_sys = sys.get_equation_systems().add_system<System>("__hypre_ads_vertices");
360 4 : lagrange_sys.hide_output() = true;
361 4 : lagrange_sys.add_variable("__lagrange");
362 4 : lagrange_sys.reinit_mesh();
363 4 : System & nedelec_sys = sys.get_equation_systems().add_system<System>("__hypre_ads_edges");
364 4 : nedelec_sys.hide_output() = true;
365 4 : nedelec_sys.add_variable("__nedelec", FIRST, NEDELEC_ONE);
366 4 : nedelec_sys.reinit_mesh();
367 :
368 : // Global (i.e. total) and local (i.e. to this processor) number of faces,
369 : // edges and vertices
370 4 : const std::vector<dof_id_type> n_all_faces = sys.get_dof_map().n_dofs_per_processor(v);
371 4 : const dof_id_type n_glb_faces = std::accumulate(n_all_faces.begin(), n_all_faces.end(), 0);
372 4 : const dof_id_type n_loc_faces = n_all_faces[global_processor_id()];
373 4 : const dof_id_type n_glb_edges = nedelec_sys.n_dofs();
374 4 : const dof_id_type n_loc_edges = nedelec_sys.n_local_dofs();
375 4 : const dof_id_type n_glb_verts = lagrange_sys.n_dofs();
376 4 : const dof_id_type n_loc_verts = lagrange_sys.n_local_dofs();
377 :
378 : // We require local indexing through the vertices and contiguous indexing
379 : // through the edges and faces: since the vertices are enumerated on their own
380 : // system, a local index can simply be obtained by subtracting those vertices
381 : // in all preceding processors; since the edges are enumerated on their own
382 : // system, we already have contiguous indexing; for the faces, if the system
383 : // we are given by the user has multiple variables/splits, we effectively need
384 : // to add local face numbers to those faces in all preceding processors
385 :
386 : // The number of vertices and faces in all preceding processors
387 0 : dof_id_signed_type vert_offset = -lagrange_sys.get_dof_map().first_dof();
388 4 : dof_id_signed_type face_offset =
389 0 : std::accumulate(n_all_faces.begin(), n_all_faces.begin() + global_processor_id(), 0);
390 :
391 : // Whether dofs are in variable-major or node-major order
392 4 : bool var_major = !libMesh::on_command_line("--node-major-dofs");
393 :
394 : // If variable-major order, we need only subtract all the preceding dofs; but
395 : // if node-major order, we need to enumerate the dofs on this processor
396 0 : std::vector<dof_id_type> idx_faces;
397 4 : if (var_major)
398 : {
399 4 : face_offset -= sys.get_dof_map().first_dof();
400 4 : for (auto i : make_range(v))
401 0 : face_offset -= sys.get_dof_map().n_local_dofs(i);
402 : }
403 : else
404 0 : sys.get_dof_map().local_variable_indices(idx_faces, sys.get_mesh(), v);
405 :
406 : // Create the discrete grandient matrix, representing the edges in terms of its vertices
407 : // Preallocate 2 diagonal + 2 off-diagonal nonzeros as the vertices could fall on either
408 4 : PetscMatrix<Real> G(Comm, n_glb_edges, n_glb_verts, n_loc_edges, n_loc_verts, 2, 2);
409 :
410 : // Create the discrete curl matrix, representing the faces in terms of its edges
411 : // Preallocate 4 diagonal + 4 off-diagonal nonzeros as the edges could fall on either
412 4 : PetscMatrix<Real> C(Comm, n_glb_faces, n_glb_edges, n_loc_faces, n_loc_edges, 4, 4);
413 :
414 : // Create array for the coordinates of the local vertices
415 : PetscReal * coords;
416 4 : LibmeshPetscCallA(comm, PetscMalloc1(dim * n_loc_verts, &coords));
417 :
418 : // Populate the discrete gradient matrix and the coordinates array
419 17130 : for (const auto & elem : sys.get_mesh().active_local_element_ptr_range())
420 49545 : for (auto face : make_range(elem->n_faces()))
421 : {
422 : // The number of edges on this face
423 40986 : const unsigned n_face_edges = Elem::type_to_n_sides_map[elem->side_type(face)];
424 :
425 : // The faces's three/four edges
426 40986 : std::vector<dof_id_type> edge_dofs(n_face_edges);
427 0 : std::vector<bool> edge_orients(n_face_edges);
428 184194 : for (auto face_edge : make_range(n_face_edges))
429 : {
430 : // Convert from face-wise to element-wise edge id
431 143208 : const unsigned edge = elem->local_side_node(face, n_face_edges + face_edge)
432 143208 : - elem->n_vertices();
433 0 : libmesh_assert(elem->is_edge_on_side(edge, face));
434 :
435 : // The edge's two vertices
436 : dof_id_type vert_dofs[2];
437 429624 : for (auto vert : make_range(2))
438 : {
439 286416 : const unsigned loc_vert_node_id = elem->local_edge_node(edge, vert);
440 0 : libmesh_assert(elem->is_vertex(loc_vert_node_id));
441 :
442 286416 : const Node & vert_node = elem->node_ref(loc_vert_node_id);
443 286416 : vert_dofs[vert] = vert_node.dof_number(lagrange_sys.number(), 0, 0);
444 :
445 : // If owned, populate coordinates array
446 286416 : if (vert_node.processor_id() == global_processor_id())
447 : {
448 0 : const dof_id_type loc_vert_dof = vert_offset + vert_dofs[vert];
449 0 : libmesh_assert(loc_vert_dof < n_loc_verts);
450 :
451 1085712 : for (auto d : make_range(dim))
452 814284 : coords[dim * loc_vert_dof + d] = vert_node(d);
453 : }
454 : }
455 :
456 : // The edge's (middle) node
457 143208 : const unsigned loc_edge_node_id = elem->local_edge_node(edge, 2);
458 0 : libmesh_assert(elem->is_edge(loc_edge_node_id));
459 :
460 143208 : const Node & edge_node = elem->node_ref(loc_edge_node_id);
461 143208 : edge_dofs[face_edge] = edge_node.dof_number(nedelec_sys.number(), 0, 0);
462 286416 : edge_orients[face_edge] = elem->positive_edge_orientation(edge) ^
463 143208 : elem->relative_edge_face_order(edge, face);
464 :
465 : // If owned, populate discrete gradient matrix
466 143208 : if (edge_node.processor_id() == global_processor_id())
467 : {
468 139630 : const Real sign = elem->positive_edge_orientation(edge) ? 1 : -1;
469 :
470 0 : libmesh_assert(edge_dofs[face_edge] >= G.row_start());
471 0 : libmesh_assert(edge_dofs[face_edge] < G.row_stop());
472 139630 : G.set(edge_dofs[face_edge], vert_dofs[0], sign);
473 139630 : G.set(edge_dofs[face_edge], vert_dofs[1], -sign);
474 : }
475 : }
476 :
477 : // The faces's (middle) node
478 40986 : const unsigned loc_face_node_id = elem->local_side_node(face, 2 * n_face_edges);
479 0 : libmesh_assert(elem->is_face(loc_face_node_id));
480 :
481 40986 : const Node & face_node = elem->node_ref(loc_face_node_id);
482 :
483 : // If owned, populate discrete curl matrix
484 40986 : if (face_node.processor_id() == global_processor_id())
485 : {
486 40580 : const dof_id_type face_dof = face_node.dof_number(sys.number(), v, 0);
487 :
488 40580 : const dof_id_type cont_face_dof = face_offset +
489 40580 : (var_major ? face_dof
490 0 : : std::distance(idx_faces.begin(),
491 : std::find(idx_faces.begin(), idx_faces.end(), face_dof)));
492 :
493 40580 : const bool face_orient = elem->positive_face_orientation(face);
494 :
495 182324 : for (auto face_edge : make_range(n_face_edges))
496 : {
497 141744 : const Real sign = face_orient ^ edge_orients[face_edge] ? 1 : -1;
498 :
499 0 : libmesh_assert(cont_face_dof >= C.row_start());
500 0 : libmesh_assert(cont_face_dof < C.row_stop());
501 141744 : C.set(cont_face_dof, edge_dofs[face_edge], sign);
502 : }
503 : }
504 : }
505 :
506 : // Assemble the discrete gradient and discrete curl matrices
507 4 : G.close();
508 4 : C.close();
509 :
510 : #ifndef NDEBUG
511 : // The product CG of the two matrices should be the zero matrix
512 0 : PetscMatrix<Real> CG(Comm);
513 0 : C.matrix_matrix_mult(G, CG);
514 0 : libmesh_assert(CG.linfty_norm() == 0);
515 : #endif
516 :
517 : // Hand over the matrices and coordinates array
518 4 : LibmeshPetscCallA(comm, PCHYPRESetDiscreteGradient(pc, G.mat()));
519 4 : LibmeshPetscCallA(comm, PCHYPRESetDiscreteCurl(pc, C.mat()));
520 4 : LibmeshPetscCallA(comm, PCSetCoordinates(pc, dim, n_loc_verts, coords));
521 :
522 : // Free allocated memory for the coordinates array
523 4 : LibmeshPetscCallA(comm, PetscFree(coords));
524 8 : }
525 : #endif
526 :
527 :
528 : //------------------------------------------------------------------
529 : // Explicit instantiations
530 : template class LIBMESH_EXPORT PetscPreconditioner<Number>;
531 :
532 : } // namespace libMesh
533 :
534 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
