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 800 : const FEType & var_type = sys.variable(v).type();
193 0 : const bool first_order_nedelec =
194 800 : var_type.order == FIRST && var_type.family == NEDELEC_ONE;
195 0 : const bool first_order_raviart_thomas =
196 800 : var_type.order == FIRST && var_type.family == RAVIART_THOMAS;
197 :
198 800 : libmesh_error_msg_if(!first_order_nedelec &&
199 : (!first_order_raviart_thomas ||
200 : sys.get_mesh().mesh_dimension() != 2),
201 : "Error applying hypre AMS to a system "
202 : "whose variable is not 1st order Nedelec or 1st "
203 : "order Raviart-Thomas on a 2d mesh");
204 800 : set_hypre_ams_data(pc, sys, v);
205 : }
206 8 : else if (std::string(hypre_type) == "ads")
207 : {
208 : // If multiple variables, we error out as senseless
209 4 : libmesh_error_msg_if(sys.n_vars() > 1,
210 : "Error applying hypre ADS to a system with multiple "
211 : "variables");
212 : // If not a 3d 1st order Raviart-Thomas system, we error out as we do not
213 : // support anything else at the moment
214 4 : const FEType & var_type = sys.variable(v).type();
215 0 : const bool first_order_raviart_thomas =
216 4 : var_type.order == FIRST && var_type.family == RAVIART_THOMAS;
217 :
218 4 : libmesh_error_msg_if(!first_order_raviart_thomas ||
219 : sys.get_mesh().mesh_dimension() != 3,
220 : "Error applying hypre ADS to a system "
221 : "whose variable is not 1st "
222 : "order Raviart-Thomas on a 3d mesh");
223 4 : set_hypre_ads_data(pc, sys, v);
224 : }
225 : }
226 259204 : }
227 :
228 :
229 :
230 : template <typename T>
231 800 : void PetscPreconditioner<T>::set_hypre_ams_data(PC & pc, System & sys, const unsigned v)
232 : {
233 : // Get the communicator from the PETSc object
234 : Parallel::communicator comm;
235 800 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
236 800 : libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
237 : "Error retrieving communicator");
238 0 : Parallel::Communicator Comm(comm);
239 :
240 : // The mesh/problem dimension
241 800 : const unsigned dim = sys.get_mesh().mesh_dimension();
242 :
243 : // Dummy Lagrange system defined over the same mesh so we can enumerate the vertices
244 800 : System & lagrange_sys = sys.get_equation_systems().add_system<System>("__hypre_ams_vertices");
245 800 : lagrange_sys.hide_output() = true;
246 800 : lagrange_sys.add_variable("__lagrange");
247 800 : lagrange_sys.reinit_mesh();
248 :
249 : // Global (i.e. total) and local (i.e. to this processor) number of edges and vertices
250 800 : const std::vector<dof_id_type> n_all_edges = sys.get_dof_map().n_dofs_per_processor(v);
251 800 : const dof_id_type n_glb_edges = std::accumulate(n_all_edges.begin(), n_all_edges.end(), 0);
252 800 : const dof_id_type n_loc_edges = n_all_edges[global_processor_id()];
253 800 : const dof_id_type n_glb_verts = lagrange_sys.n_dofs();
254 800 : const dof_id_type n_loc_verts = lagrange_sys.n_local_dofs();
255 :
256 : // We require local indexing through the vertices and contiguous indexing
257 : // through the edges: since the vertices are enumerated on their own system,
258 : // a local index can simply be obtained by subtracting those vertices in all
259 : // preceding processors; for the edges, if the system we are given by the
260 : // user has multiple variables/splits, we effectively need to add local edge
261 : // numbers to those edges in all preceding processors
262 :
263 : // The number of vertices and edges in all preceding processors
264 0 : dof_id_signed_type vert_offset = -lagrange_sys.get_dof_map().first_dof();
265 800 : dof_id_signed_type edge_offset =
266 0 : std::accumulate(n_all_edges.begin(), n_all_edges.begin() + global_processor_id(), 0);
267 :
268 : // Whether dofs are in variable-major or node-major order
269 800 : bool var_major = !libMesh::on_command_line("--node-major-dofs");
270 :
271 : // If variable-major order, we need only subtract all the preceding dofs; but
272 : // if node-major order, we need to enumerate the dofs on this processor
273 0 : std::vector<dof_id_type> idx_edges;
274 800 : if (var_major)
275 : {
276 800 : edge_offset -= sys.get_dof_map().first_dof();
277 800 : for (auto i : make_range(v))
278 0 : edge_offset -= sys.get_dof_map().n_local_dofs(i);
279 : }
280 : else
281 0 : sys.get_dof_map().local_variable_indices(idx_edges, sys.get_mesh(), v);
282 :
283 : // Create the discrete grandient matrix, representing the edges in terms of its vertices
284 : // Preallocate 2 diagonal + 2 off-diagonal nonzeros as the vertices could fall on either
285 800 : PetscMatrix<Real> G(Comm, n_glb_edges, n_glb_verts, n_loc_edges, n_loc_verts, 2, 2);
286 :
287 : // Create array for the coordinates of the local vertices
288 : PetscReal * coords;
289 800 : LibmeshPetscCallA(comm, PetscMalloc1(dim * n_loc_verts, &coords));
290 :
291 : // Populate the discrete gradient matrix and the coordinates array
292 149100 : for (const auto & elem : sys.get_mesh().active_local_element_ptr_range())
293 401562 : for (auto edge : make_range(elem->n_edges()))
294 : {
295 : // The edge's two vertices
296 : dof_id_type vert_dofs[2];
297 984636 : for (auto vert : make_range(2))
298 : {
299 656424 : const unsigned loc_vert_node_id = elem->local_edge_node(edge, vert);
300 0 : libmesh_assert(elem->is_vertex(loc_vert_node_id));
301 :
302 656424 : const Node & vert_node = elem->node_ref(loc_vert_node_id);
303 656424 : vert_dofs[vert] = vert_node.dof_number(lagrange_sys.number(), 0, 0);
304 :
305 : // If owned, populate coordinates array
306 656424 : if (vert_node.processor_id() == global_processor_id())
307 : {
308 0 : const dof_id_type loc_vert_dof = vert_offset + vert_dofs[vert];
309 0 : libmesh_assert(loc_vert_dof < n_loc_verts);
310 :
311 1902888 : for (auto d : make_range(dim))
312 1360376 : coords[dim * loc_vert_dof + d] = vert_node(d);
313 : }
314 : }
315 :
316 : // The edge's (middle) node
317 328212 : const unsigned loc_edge_node_id = elem->local_edge_node(edge, 2);
318 0 : libmesh_assert(elem->is_edge(loc_edge_node_id));
319 :
320 328212 : const Node & edge_node = elem->node_ref(loc_edge_node_id);
321 :
322 : // If owned, populate discrete gradient matrix
323 328212 : if (edge_node.processor_id() == global_processor_id())
324 : {
325 303502 : const dof_id_type edge_dof = edge_node.dof_number(sys.number(), v, 0);
326 :
327 303502 : const dof_id_type cont_edge_dof = edge_offset +
328 303502 : (var_major ? edge_dof
329 0 : : std::distance(idx_edges.begin(),
330 : std::find(idx_edges.begin(), idx_edges.end(), edge_dof)));
331 :
332 303502 : const Real sign = elem->positive_edge_orientation(edge) ? 1 : -1;
333 :
334 0 : libmesh_assert(cont_edge_dof >= G.row_start());
335 0 : libmesh_assert(cont_edge_dof < G.row_stop());
336 303502 : G.set(cont_edge_dof, vert_dofs[0], sign);
337 303502 : G.set(cont_edge_dof, vert_dofs[1], -sign);
338 : }
339 : }
340 :
341 : // Assemble the discrete gradient matrix
342 800 : G.close();
343 :
344 : // Hand over the matrix and coordinates array
345 800 : LibmeshPetscCallA(comm, PCHYPRESetDiscreteGradient(pc, G.mat()));
346 800 : LibmeshPetscCallA(comm, PCSetCoordinates(pc, dim, n_loc_verts, coords));
347 :
348 : // Free allocated memory for the coordinates array
349 800 : LibmeshPetscCallA(comm, PetscFree(coords));
350 1600 : }
351 :
352 :
353 :
354 : template <typename T>
355 4 : void PetscPreconditioner<T>::set_hypre_ads_data(PC & pc, System & sys, const unsigned v)
356 : {
357 : // Get the communicator from the PETSc object
358 : Parallel::communicator comm;
359 4 : PetscErrorCode ierr = PetscObjectGetComm((PetscObject)pc, &comm);
360 4 : libmesh_error_msg_if(ierr != LIBMESH_PETSC_SUCCESS,
361 : "Error retrieving communicator");
362 0 : Parallel::Communicator Comm(comm);
363 :
364 : // The mesh/problem dimension
365 4 : const unsigned dim = sys.get_mesh().mesh_dimension();
366 :
367 : // Dummy Lagrange and Nédélec systems defined over the same mesh so we can
368 : // enumerate the vertices and edges, respectively
369 4 : System & lagrange_sys = sys.get_equation_systems().add_system<System>("__hypre_ads_vertices");
370 4 : lagrange_sys.hide_output() = true;
371 4 : lagrange_sys.add_variable("__lagrange");
372 4 : lagrange_sys.reinit_mesh();
373 4 : System & nedelec_sys = sys.get_equation_systems().add_system<System>("__hypre_ads_edges");
374 4 : nedelec_sys.hide_output() = true;
375 4 : nedelec_sys.add_variable("__nedelec", FIRST, NEDELEC_ONE);
376 4 : nedelec_sys.reinit_mesh();
377 :
378 : // Global (i.e. total) and local (i.e. to this processor) number of faces,
379 : // edges and vertices
380 4 : const std::vector<dof_id_type> n_all_faces = sys.get_dof_map().n_dofs_per_processor(v);
381 4 : const dof_id_type n_glb_faces = std::accumulate(n_all_faces.begin(), n_all_faces.end(), 0);
382 4 : const dof_id_type n_loc_faces = n_all_faces[global_processor_id()];
383 4 : const dof_id_type n_glb_edges = nedelec_sys.n_dofs();
384 4 : const dof_id_type n_loc_edges = nedelec_sys.n_local_dofs();
385 4 : const dof_id_type n_glb_verts = lagrange_sys.n_dofs();
386 4 : const dof_id_type n_loc_verts = lagrange_sys.n_local_dofs();
387 :
388 : // We require local indexing through the vertices and contiguous indexing
389 : // through the edges and faces: since the vertices are enumerated on their own
390 : // system, a local index can simply be obtained by subtracting those vertices
391 : // in all preceding processors; since the edges are enumerated on their own
392 : // system, we already have contiguous indexing; for the faces, if the system
393 : // we are given by the user has multiple variables/splits, we effectively need
394 : // to add local face numbers to those faces in all preceding processors
395 :
396 : // The number of vertices and faces in all preceding processors
397 0 : dof_id_signed_type vert_offset = -lagrange_sys.get_dof_map().first_dof();
398 4 : dof_id_signed_type face_offset =
399 0 : std::accumulate(n_all_faces.begin(), n_all_faces.begin() + global_processor_id(), 0);
400 :
401 : // Whether dofs are in variable-major or node-major order
402 4 : bool var_major = !libMesh::on_command_line("--node-major-dofs");
403 :
404 : // If variable-major order, we need only subtract all the preceding dofs; but
405 : // if node-major order, we need to enumerate the dofs on this processor
406 0 : std::vector<dof_id_type> idx_faces;
407 4 : if (var_major)
408 : {
409 4 : face_offset -= sys.get_dof_map().first_dof();
410 4 : for (auto i : make_range(v))
411 0 : face_offset -= sys.get_dof_map().n_local_dofs(i);
412 : }
413 : else
414 0 : sys.get_dof_map().local_variable_indices(idx_faces, sys.get_mesh(), v);
415 :
416 : // Create the discrete grandient matrix, representing the edges in terms of its vertices
417 : // Preallocate 2 diagonal + 2 off-diagonal nonzeros as the vertices could fall on either
418 4 : PetscMatrix<Real> G(Comm, n_glb_edges, n_glb_verts, n_loc_edges, n_loc_verts, 2, 2);
419 :
420 : // Create the discrete curl matrix, representing the faces in terms of its edges
421 : // Preallocate 4 diagonal + 4 off-diagonal nonzeros as the edges could fall on either
422 4 : PetscMatrix<Real> C(Comm, n_glb_faces, n_glb_edges, n_loc_faces, n_loc_edges, 4, 4);
423 :
424 : // Create array for the coordinates of the local vertices
425 : PetscReal * coords;
426 4 : LibmeshPetscCallA(comm, PetscMalloc1(dim * n_loc_verts, &coords));
427 :
428 : // Populate the discrete gradient matrix and the coordinates array
429 17130 : for (const auto & elem : sys.get_mesh().active_local_element_ptr_range())
430 49545 : for (auto face : make_range(elem->n_faces()))
431 : {
432 : // The number of edges on this face
433 40986 : const unsigned n_face_edges = Elem::type_to_n_sides_map[elem->side_type(face)];
434 :
435 : // The faces's three/four edges
436 40986 : std::vector<dof_id_type> edge_dofs(n_face_edges);
437 0 : std::vector<bool> edge_orients(n_face_edges);
438 184194 : for (auto face_edge : make_range(n_face_edges))
439 : {
440 : // Convert from face-wise to element-wise edge id
441 143208 : const unsigned edge = elem->local_side_node(face, n_face_edges + face_edge)
442 143208 : - elem->n_vertices();
443 0 : libmesh_assert(elem->is_edge_on_side(edge, face));
444 :
445 : // The edge's two vertices
446 : dof_id_type vert_dofs[2];
447 429624 : for (auto vert : make_range(2))
448 : {
449 286416 : const unsigned loc_vert_node_id = elem->local_edge_node(edge, vert);
450 0 : libmesh_assert(elem->is_vertex(loc_vert_node_id));
451 :
452 286416 : const Node & vert_node = elem->node_ref(loc_vert_node_id);
453 286416 : vert_dofs[vert] = vert_node.dof_number(lagrange_sys.number(), 0, 0);
454 :
455 : // If owned, populate coordinates array
456 286416 : if (vert_node.processor_id() == global_processor_id())
457 : {
458 0 : const dof_id_type loc_vert_dof = vert_offset + vert_dofs[vert];
459 0 : libmesh_assert(loc_vert_dof < n_loc_verts);
460 :
461 1085712 : for (auto d : make_range(dim))
462 814284 : coords[dim * loc_vert_dof + d] = vert_node(d);
463 : }
464 : }
465 :
466 : // The edge's (middle) node
467 143208 : const unsigned loc_edge_node_id = elem->local_edge_node(edge, 2);
468 0 : libmesh_assert(elem->is_edge(loc_edge_node_id));
469 :
470 143208 : const Node & edge_node = elem->node_ref(loc_edge_node_id);
471 143208 : edge_dofs[face_edge] = edge_node.dof_number(nedelec_sys.number(), 0, 0);
472 286416 : edge_orients[face_edge] = elem->positive_edge_orientation(edge) ^
473 143208 : elem->relative_edge_face_order(edge, face);
474 :
475 : // If owned, populate discrete gradient matrix
476 143208 : if (edge_node.processor_id() == global_processor_id())
477 : {
478 139630 : const Real sign = elem->positive_edge_orientation(edge) ? 1 : -1;
479 :
480 0 : libmesh_assert(edge_dofs[face_edge] >= G.row_start());
481 0 : libmesh_assert(edge_dofs[face_edge] < G.row_stop());
482 139630 : G.set(edge_dofs[face_edge], vert_dofs[0], sign);
483 139630 : G.set(edge_dofs[face_edge], vert_dofs[1], -sign);
484 : }
485 : }
486 :
487 : // The faces's (middle) node
488 40986 : const unsigned loc_face_node_id = elem->local_side_node(face, 2 * n_face_edges);
489 0 : libmesh_assert(elem->is_face(loc_face_node_id));
490 :
491 40986 : const Node & face_node = elem->node_ref(loc_face_node_id);
492 :
493 : // If owned, populate discrete curl matrix
494 40986 : if (face_node.processor_id() == global_processor_id())
495 : {
496 40580 : const dof_id_type face_dof = face_node.dof_number(sys.number(), v, 0);
497 :
498 40580 : const dof_id_type cont_face_dof = face_offset +
499 40580 : (var_major ? face_dof
500 0 : : std::distance(idx_faces.begin(),
501 : std::find(idx_faces.begin(), idx_faces.end(), face_dof)));
502 :
503 40580 : const bool face_orient = elem->positive_face_orientation(face);
504 :
505 182324 : for (auto face_edge : make_range(n_face_edges))
506 : {
507 141744 : const Real sign = face_orient ^ edge_orients[face_edge] ? 1 : -1;
508 :
509 0 : libmesh_assert(cont_face_dof >= C.row_start());
510 0 : libmesh_assert(cont_face_dof < C.row_stop());
511 141744 : C.set(cont_face_dof, edge_dofs[face_edge], sign);
512 : }
513 : }
514 : }
515 :
516 : // Assemble the discrete gradient and discrete curl matrices
517 4 : G.close();
518 4 : C.close();
519 :
520 : #ifndef NDEBUG
521 : // The product CG of the two matrices should be the zero matrix
522 0 : PetscMatrix<Real> CG(Comm);
523 0 : C.matrix_matrix_mult(G, CG);
524 0 : libmesh_assert(CG.linfty_norm() == 0);
525 : #endif
526 :
527 : // Hand over the matrices and coordinates array
528 4 : LibmeshPetscCallA(comm, PCHYPRESetDiscreteGradient(pc, G.mat()));
529 4 : LibmeshPetscCallA(comm, PCHYPRESetDiscreteCurl(pc, C.mat()));
530 4 : LibmeshPetscCallA(comm, PCSetCoordinates(pc, dim, n_loc_verts, coords));
531 :
532 : // Free allocated memory for the coordinates array
533 4 : LibmeshPetscCallA(comm, PetscFree(coords));
534 8 : }
535 : #endif
536 :
537 :
538 : //------------------------------------------------------------------
539 : // Explicit instantiations
540 : template class LIBMESH_EXPORT PetscPreconditioner<Number>;
541 :
542 : } // namespace libMesh
543 :
544 : #endif // #ifdef LIBMESH_HAVE_PETSC
|
