MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
bilininteg_diffusion_patch.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
11
12
13#include "../fem.hpp"
14#include "../../mesh/nurbs.hpp"
15
16#include "../../linalg/dtensor.hpp" // For Reshape
17#include "../../general/forall.hpp"
18
19using namespace std;
20
21namespace mfem
22{
23
24// Adapted from PADiffusionSetup3D
25void SetupPatch3D(const int Q1Dx,
26 const int Q1Dy,
27 const int Q1Dz,
28 const int coeffDim,
29 const bool symmetric,
30 const Array<real_t> &w,
31 const Vector &j,
32 const Vector &c,
33 Vector &d)
34{
35 const bool const_c = (c.Size() == 1);
36 MFEM_VERIFY(coeffDim < 6 ||
37 !const_c, "Constant matrix coefficient not supported");
38
39 const auto W = Reshape(w.Read(), Q1Dx,Q1Dy,Q1Dz);
40 const auto J = Reshape(j.Read(), Q1Dx,Q1Dy,Q1Dz,3,3);
41 const auto C = const_c ? Reshape(c.Read(), 1,1,1,1) :
42 Reshape(c.Read(), coeffDim,Q1Dx,Q1Dy,Q1Dz);
43 d.SetSize(Q1Dx * Q1Dy * Q1Dz * (symmetric ? 6 : 9));
44 auto D = Reshape(d.Write(), Q1Dx,Q1Dy,Q1Dz, symmetric ? 6 : 9);
45 const int NE = 1; // TODO: MFEM_FORALL_3D without e?
46 MFEM_FORALL_3D(e, NE, Q1Dx, Q1Dy, Q1Dz,
47 {
48 MFEM_FOREACH_THREAD(qx,x,Q1Dx)
49 {
50 MFEM_FOREACH_THREAD(qy,y,Q1Dy)
51 {
52 MFEM_FOREACH_THREAD(qz,z,Q1Dz)
53 {
54 const real_t J11 = J(qx,qy,qz,0,0);
55 const real_t J21 = J(qx,qy,qz,1,0);
56 const real_t J31 = J(qx,qy,qz,2,0);
57 const real_t J12 = J(qx,qy,qz,0,1);
58 const real_t J22 = J(qx,qy,qz,1,1);
59 const real_t J32 = J(qx,qy,qz,2,1);
60 const real_t J13 = J(qx,qy,qz,0,2);
61 const real_t J23 = J(qx,qy,qz,1,2);
62 const real_t J33 = J(qx,qy,qz,2,2);
63 const real_t detJ = J11 * (J22 * J33 - J32 * J23) -
64 /* */ J21 * (J12 * J33 - J32 * J13) +
65 /* */ J31 * (J12 * J23 - J22 * J13);
66 const real_t w_detJ = W(qx,qy,qz) / detJ;
67 // adj(J)
68 const real_t A11 = (J22 * J33) - (J23 * J32);
69 const real_t A12 = (J32 * J13) - (J12 * J33);
70 const real_t A13 = (J12 * J23) - (J22 * J13);
71 const real_t A21 = (J31 * J23) - (J21 * J33);
72 const real_t A22 = (J11 * J33) - (J13 * J31);
73 const real_t A23 = (J21 * J13) - (J11 * J23);
74 const real_t A31 = (J21 * J32) - (J31 * J22);
75 const real_t A32 = (J31 * J12) - (J11 * J32);
76 const real_t A33 = (J11 * J22) - (J12 * J21);
77
78 if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
79 {
80 // Compute entries of R = MJ^{-T} = M adj(J)^T, without det J.
81 const real_t M11 = C(0, qx,qy,qz);
82 const real_t M12 = C(1, qx,qy,qz);
83 const real_t M13 = C(2, qx,qy,qz);
84 const real_t M21 = (!symmetric) ? C(3, qx,qy,qz) : M12;
85 const real_t M22 = (!symmetric) ? C(4, qx,qy,qz) : C(3, qx,qy,qz);
86 const real_t M23 = (!symmetric) ? C(5, qx,qy,qz) : C(4, qx,qy,qz);
87 const real_t M31 = (!symmetric) ? C(6, qx,qy,qz) : M13;
88 const real_t M32 = (!symmetric) ? C(7, qx,qy,qz) : M23;
89 const real_t M33 = (!symmetric) ? C(8, qx,qy,qz) : C(5, qx,qy,qz);
90
91 const real_t R11 = M11*A11 + M12*A12 + M13*A13;
92 const real_t R12 = M11*A21 + M12*A22 + M13*A23;
93 const real_t R13 = M11*A31 + M12*A32 + M13*A33;
94 const real_t R21 = M21*A11 + M22*A12 + M23*A13;
95 const real_t R22 = M21*A21 + M22*A22 + M23*A23;
96 const real_t R23 = M21*A31 + M22*A32 + M23*A33;
97 const real_t R31 = M31*A11 + M32*A12 + M33*A13;
98 const real_t R32 = M31*A21 + M32*A22 + M33*A23;
99 const real_t R33 = M31*A31 + M32*A32 + M33*A33;
100
101 // Now set D to J^{-1} R = adj(J) R
102 D(qx,qy,qz,0) = w_detJ * (A11*R11 + A12*R21 + A13*R31); // 1,1
103 const real_t D12 = w_detJ * (A11*R12 + A12*R22 + A13*R32);
104 D(qx,qy,qz,1) = D12; // 1,2
105 D(qx,qy,qz,2) = w_detJ * (A11*R13 + A12*R23 + A13*R33); // 1,3
106
107 const real_t D21 = w_detJ * (A21*R11 + A22*R21 + A23*R31);
108 const real_t D22 = w_detJ * (A21*R12 + A22*R22 + A23*R32);
109 const real_t D23 = w_detJ * (A21*R13 + A22*R23 + A23*R33);
110
111 const real_t D33 = w_detJ * (A31*R13 + A32*R23 + A33*R33);
112
113 D(qx,qy,qz,3) = symmetric ? D22 : D21; // 2,2 or 2,1
114 D(qx,qy,qz,4) = symmetric ? D23 : D22; // 2,3 or 2,2
115 D(qx,qy,qz,5) = symmetric ? D33 : D23; // 3,3 or 2,3
116
117 if (!symmetric)
118 {
119 D(qx,qy,qz,6) = w_detJ * (A31*R11 + A32*R21 + A33*R31); // 3,1
120 D(qx,qy,qz,7) = w_detJ * (A31*R12 + A32*R22 + A33*R32); // 3,2
121 D(qx,qy,qz,8) = D33; // 3,3
122 }
123 }
124 else // Vector or scalar coefficient version
125 {
126 const real_t C1 = const_c ? C(0,0,0,0) : C(0,qx,qy,qz);
127 const real_t C2 = const_c ? C(0,0,0,0) :
128 (coeffDim == 3 ? C(1,qx,qy,qz) : C(0,qx,qy,qz));
129 const real_t C3 = const_c ? C(0,0,0,0) :
130 (coeffDim == 3 ? C(2,qx,qy,qz) : C(0,qx,qy,qz));
131
132 // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
133 D(qx,qy,qz,0) = w_detJ * (C1*A11*A11 + C2*A12*A12 + C3*A13*A13); // 1,1
134 D(qx,qy,qz,1) = w_detJ * (C1*A11*A21 + C2*A12*A22 + C3*A13*A23); // 2,1
135 D(qx,qy,qz,2) = w_detJ * (C1*A11*A31 + C2*A12*A32 + C3*A13*A33); // 3,1
136 D(qx,qy,qz,3) = w_detJ * (C1*A21*A21 + C2*A22*A22 + C3*A23*A23); // 2,2
137 D(qx,qy,qz,4) = w_detJ * (C1*A21*A31 + C2*A22*A32 + C3*A23*A33); // 3,2
138 D(qx,qy,qz,5) = w_detJ * (C1*A31*A31 + C2*A32*A32 + C3*A33*A33); // 3,3
139 }
140 }
141 }
142 }
143 });
144}
145
146// Compute a reduced integration rule, using NNLSSolver, for DiffusionIntegrator
147// on a NURBS patch with partial assembly.
148void GetReducedRule(const int nq, const int nd,
149 Array2D<real_t> const& B,
150 Array2D<real_t> const& G,
151 std::vector<int> minQ,
152 std::vector<int> maxQ,
153 std::vector<int> minD,
154 std::vector<int> maxD,
155 std::vector<int> minDD,
156 std::vector<int> maxDD,
157 const IntegrationRule *ir,
158 const bool zeroOrder,
159 std::vector<Vector> & reducedWeights,
160 std::vector<std::vector<int>> & reducedIDs)
161{
162 MFEM_VERIFY(B.NumRows() == nq, "");
163 MFEM_VERIFY(B.NumCols() == nd, "");
164 MFEM_VERIFY(G.NumRows() == nq, "");
165 MFEM_VERIFY(G.NumCols() == nd, "");
166 MFEM_VERIFY(ir->GetNPoints() == nq, "");
167
168 for (int dof=0; dof<nd; ++dof)
169 {
170 // Integrate diffusion for B(:,dof) against all other B(:,i)
171
172 const int nc_dof = maxDD[dof] - minDD[dof] + 1;
173 const int nw_dof = maxD[dof] - minD[dof] + 1;
174
175 // G is of size nc_dof x nw_dof
176 MFEM_VERIFY(nc_dof <= nw_dof, "The NNLS system for the reduced "
177 "integration rule requires more full integration points. Try"
178 " increasing the order of the full integration rule.");
179 DenseMatrix Gmat(nc_dof, nw_dof);
180 Gmat = 0.0;
181
182 Vector w(nw_dof);
183 w = 0.0;
184
185 for (int qx = minD[dof]; qx <= maxD[dof]; ++qx)
186 {
187 const real_t Bq = zeroOrder ? B(qx,dof) : G(qx,dof);
188
189 const IntegrationPoint &ip = ir->IntPoint(qx);
190 const real_t w_qx = ip.weight;
191 w[qx - minD[dof]] = w_qx;
192
193 for (int dx = minQ[qx]; dx <= maxQ[qx]; ++dx)
194 {
195 const real_t Bd = zeroOrder ? B(qx,dx) : G(qx,dx);
196
197 Gmat(dx - minDD[dof], qx - minD[dof]) = Bq * Bd;
198 }
199 }
200
201 Vector sol(Gmat.NumCols());
202
203#ifdef MFEM_USE_LAPACK
204 NNLSSolver nnls;
205 nnls.SetOperator(Gmat);
206
207 nnls.Mult(w, sol);
208#else
209 MFEM_ABORT("NNLSSolver requires building with LAPACK");
210#endif
211
212 int nnz = 0;
213 for (int i=0; i<sol.Size(); ++i)
214 {
215 if (sol(i) != 0.0)
216 {
217 nnz++;
218 }
219 }
220
221 MFEM_VERIFY(nnz > 0, "");
222
223 Vector wred(nnz);
224 std::vector<int> idnnz(nnz);
225 nnz = 0;
226 for (int i=0; i<sol.Size(); ++i)
227 {
228 if (sol(i) != 0.0)
229 {
230 wred[nnz] = sol[i];
231 idnnz[nnz] = i;
232 nnz++;
233 }
234 }
235
236 reducedWeights.push_back(wred);
237 reducedIDs.push_back(idnnz);
238 }
239}
240
241// Adapted from AssemblePA
242void DiffusionIntegrator::SetupPatchPA(const int patch, Mesh *mesh,
243 bool unitWeights)
244{
245 const Array<int>& Q1D = pQ1D[patch];
246 const Array<int>& D1D = pD1D[patch];
247 const std::vector<Array2D<real_t>>& B = pB[patch];
248 const std::vector<Array2D<real_t>>& G = pG[patch];
249
250 const IntArrayVar2D& minD = pminD[patch];
251 const IntArrayVar2D& maxD = pmaxD[patch];
252 const IntArrayVar2D& minQ = pminQ[patch];
253 const IntArrayVar2D& maxQ = pmaxQ[patch];
254
255 const IntArrayVar2D& minDD = pminDD[patch];
256 const IntArrayVar2D& maxDD = pmaxDD[patch];
257
258 const Array<const IntegrationRule*>& ir1d = pir1d[patch];
259
260 MFEM_VERIFY(Q1D.Size() == 3, "");
261
262 const int dims = dim; // TODO: generalize
263 const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
264
265 int nq = Q1D[0];
266 for (int i=1; i<dim; ++i)
267 {
268 nq *= Q1D[i];
269 }
270
271 int coeffDim = 1;
272 Vector coeff;
273 Array<real_t> weights(nq);
274 const int MQfullDim = MQ ? MQ->GetHeight() * MQ->GetWidth() : 0;
275 IntegrationPoint ip;
276
277 Vector jac(dim * dim * nq); // Computed as in GeometricFactors::Compute
278
279 for (int qz=0; qz<Q1D[2]; ++qz)
280 {
281 for (int qy=0; qy<Q1D[1]; ++qy)
282 {
283 for (int qx=0; qx<Q1D[0]; ++qx)
284 {
285 const int p = qx + (qy * Q1D[0]) + (qz * Q1D[0] * Q1D[1]);
286 patchRules->GetIntegrationPointFrom1D(patch, qx, qy, qz, ip);
287 const int e = patchRules->GetPointElement(patch, qx, qy, qz);
288 ElementTransformation *tr = mesh->GetElementTransformation(e);
289
290 weights[p] = ip.weight;
291
292 tr->SetIntPoint(&ip);
293
294 const DenseMatrix& Jp = tr->Jacobian();
295 for (int i=0; i<dim; ++i)
296 for (int j=0; j<dim; ++j)
297 {
298 jac[p + ((i + (j * dim)) * nq)] = Jp(i,j);
299 }
300 }
301 }
302 }
303
304 if (auto *SMQ = dynamic_cast<SymmetricMatrixCoefficient *>(MQ))
305 {
306 MFEM_VERIFY(SMQ->GetSize() == dim, "");
307 coeffDim = symmDims;
308 coeff.SetSize(symmDims * nq);
309
310 DenseSymmetricMatrix sym_mat;
311 sym_mat.SetSize(dim);
312
313 auto C = Reshape(coeff.HostWrite(), symmDims, nq);
314 for (int qz=0; qz<Q1D[2]; ++qz)
315 {
316 for (int qy=0; qy<Q1D[1]; ++qy)
317 {
318 for (int qx=0; qx<Q1D[0]; ++qx)
319 {
320 const int p = qx + (qy * Q1D[0]) + (qz * Q1D[0] * Q1D[1]);
321 const int e = patchRules->GetPointElement(patch, qx, qy, qz);
322 ElementTransformation *tr = mesh->GetElementTransformation(e);
323 patchRules->GetIntegrationPointFrom1D(patch, qx, qy, qz, ip);
324
325 SMQ->Eval(sym_mat, *tr, ip);
326 int cnt = 0;
327 for (int i=0; i<dim; ++i)
328 for (int j=i; j<dim; ++j, ++cnt)
329 {
330 C(cnt, p) = sym_mat(i,j);
331 }
332 }
333 }
334 }
335 }
336 else if (MQ)
337 {
338 symmetric = false;
339 MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
340
341 coeffDim = MQfullDim;
342
343 coeff.SetSize(MQfullDim * nq);
344
345 DenseMatrix mat;
346 mat.SetSize(dim);
347
348 auto C = Reshape(coeff.HostWrite(), MQfullDim, nq);
349 for (int qz=0; qz<Q1D[2]; ++qz)
350 {
351 for (int qy=0; qy<Q1D[1]; ++qy)
352 {
353 for (int qx=0; qx<Q1D[0]; ++qx)
354 {
355 const int p = qx + (qy * Q1D[0]) + (qz * Q1D[0] * Q1D[1]);
356 const int e = patchRules->GetPointElement(patch, qx, qy, qz);
357 ElementTransformation *tr = mesh->GetElementTransformation(e);
358 patchRules->GetIntegrationPointFrom1D(patch, qx, qy, qz, ip);
359
360 MQ->Eval(mat, *tr, ip);
361 for (int i=0; i<dim; ++i)
362 for (int j=0; j<dim; ++j)
363 {
364 C(j+(i*dim), p) = mat(i,j);
365 }
366 }
367 }
368 }
369 }
370 else if (VQ)
371 {
372 MFEM_VERIFY(VQ->GetVDim() == dim, "");
373 coeffDim = VQ->GetVDim();
374 coeff.SetSize(coeffDim * nq);
375 auto C = Reshape(coeff.HostWrite(), coeffDim, nq);
376 Vector DM(coeffDim);
377 for (int qz=0; qz<Q1D[2]; ++qz)
378 {
379 for (int qy=0; qy<Q1D[1]; ++qy)
380 {
381 for (int qx=0; qx<Q1D[0]; ++qx)
382 {
383 const int p = qx + (qy * Q1D[0]) + (qz * Q1D[0] * Q1D[1]);
384 const int e = patchRules->GetPointElement(patch, qx, qy, qz);
385 ElementTransformation *tr = mesh->GetElementTransformation(e);
386 patchRules->GetIntegrationPointFrom1D(patch, qx, qy, qz, ip);
387
388 VQ->Eval(DM, *tr, ip);
389 for (int i=0; i<coeffDim; ++i)
390 {
391 C(i, p) = DM[i];
392 }
393 }
394 }
395 }
396 }
397 else if (Q == nullptr)
398 {
399 coeff.SetSize(1);
400 coeff(0) = 1.0;
401 }
402 else if (ConstantCoefficient* cQ = dynamic_cast<ConstantCoefficient*>(Q))
403 {
404 coeff.SetSize(1);
405 coeff(0) = cQ->constant;
406 }
407 else if (dynamic_cast<QuadratureFunctionCoefficient*>(Q))
408 {
409 MFEM_ABORT("QuadratureFunction not supported yet\n");
410 }
411 else
412 {
413 coeff.SetSize(nq);
414 auto C = Reshape(coeff.HostWrite(), nq);
415 for (int qz=0; qz<Q1D[2]; ++qz)
416 {
417 for (int qy=0; qy<Q1D[1]; ++qy)
418 {
419 for (int qx=0; qx<Q1D[0]; ++qx)
420 {
421 const int p = qx + (qy * Q1D[0]) + (qz * Q1D[0] * Q1D[1]);
422 const int e = patchRules->GetPointElement(patch, qx, qy, qz);
423 ElementTransformation *tr = mesh->GetElementTransformation(e);
424 patchRules->GetIntegrationPointFrom1D(patch, qx, qy, qz, ip);
425
426 C(p) = Q->Eval(*tr, ip);
427 }
428 }
429 }
430 }
431
432 if (unitWeights)
433 {
434 weights = 1.0;
435 }
436
437 SetupPatch3D(Q1D[0], Q1D[1], Q1D[2], coeffDim, symmetric, weights, jac,
438 coeff, pa_data);
439
440 numPatches = mesh->NURBSext->GetNP();
441
442 if (integrationMode != PATCHWISE_REDUCED)
443 {
444 return;
445 }
446
447 // Solve for reduced 1D quadrature rules
448 const int totalDim = numPatches * dim * numTypes;
449 reducedWeights.resize(totalDim);
450 reducedIDs.resize(totalDim);
451
452 auto rw = Reshape(reducedWeights.data(), numTypes, dim, numPatches);
453 auto rid = Reshape(reducedIDs.data(), numTypes, dim, numPatches);
454
455 for (int d=0; d<dim; ++d)
456 {
457 // The reduced rules could be cached to avoid repeated computation, but
458 // the cost of this setup seems low.
459 GetReducedRule(Q1D[d], D1D[d], B[d], G[d],
460 minQ[d], maxQ[d],
461 minD[d], maxD[d],
462 minDD[d], maxDD[d], ir1d[d], true,
463 rw(0,d,patch), rid(0,d,patch));
464 GetReducedRule(Q1D[d], D1D[d], B[d], G[d],
465 minQ[d], maxQ[d],
466 minD[d], maxD[d],
467 minDD[d], maxDD[d], ir1d[d], false,
468 rw(1,d,patch), rid(1,d,patch));
469 }
470}
471
472// This version uses full 1D quadrature rules, taking into account the
473// minimum interaction between basis functions and integration points.
474void DiffusionIntegrator::AssemblePatchMatrix_fullQuadrature(
475 const int patch, const FiniteElementSpace &fes, SparseMatrix*& smat)
476{
477 MFEM_VERIFY(patchRules, "patchRules must be defined");
478 dim = patchRules->GetDim();
479 const int spaceDim = dim; // TODO: generalize?
480
481 Mesh *mesh = fes.GetMesh();
482
483 if (VQ)
484 {
485 MFEM_VERIFY(VQ->GetVDim() == spaceDim,
486 "Unexpected dimension for VectorCoefficient");
487 }
488 if (MQ)
489 {
490 MFEM_VERIFY(MQ->GetWidth() == spaceDim,
491 "Unexpected width for MatrixCoefficient");
492 MFEM_VERIFY(MQ->GetHeight() == spaceDim,
493 "Unexpected height for MatrixCoefficient");
494 }
495
496#ifdef MFEM_THREAD_SAFE
497 DenseMatrix M(MQ ? spaceDim : 0);
498 Vector D(VQ ? VQ->GetVDim() : 0);
499#else
500 M.SetSize(MQ ? spaceDim : 0);
501 D.SetSize(VQ ? VQ->GetVDim() : 0);
502#endif
503
504 SetupPatchBasisData(mesh, patch);
505
506 SetupPatchPA(patch, mesh);
507
508 const Array<int>& Q1D = pQ1D[patch];
509 const Array<int>& D1D = pD1D[patch];
510 const std::vector<Array2D<real_t>>& B = pB[patch];
511 const std::vector<Array2D<real_t>>& G = pG[patch];
512
513 const IntArrayVar2D& minD = pminD[patch];
514 const IntArrayVar2D& maxD = pmaxD[patch];
515 const IntArrayVar2D& minQ = pminQ[patch];
516 const IntArrayVar2D& maxQ = pmaxQ[patch];
517
518 const IntArrayVar2D& minDD = pminDD[patch];
519 const IntArrayVar2D& maxDD = pmaxDD[patch];
520
521 int ndof = D1D[0];
522 for (int d=1; d<dim; ++d)
523 {
524 ndof *= D1D[d];
525 }
526
527 MFEM_VERIFY(3 == dim, "Only 3D so far");
528
529 // Setup quadrature point data.
530 const auto qd = Reshape(pa_data.Read(), Q1D[0]*Q1D[1]*Q1D[2],
531 (symmetric ? 6 : 9));
532
533 // NOTE: the following is adapted from PADiffusionApply3D.
534 std::vector<Array3D<real_t>> grad(dim);
535
536 for (int d=0; d<dim; ++d)
537 {
538 grad[d].SetSize(Q1D[0], Q1D[1], Q1D[2]);
539 }
540
541 Array3D<real_t> gradDXY(D1D[0], D1D[1], dim);
542 Array2D<real_t> gradDX(D1D[0], dim);
543
544 int nd[3];
545 Array3D<int> cdofs;
546
547 int *smati = nullptr;
548 int *smatj = nullptr;
549 real_t *smata = nullptr;
550 int nnz = 0;
551
552 Array<int> maxw(dim);
553 maxw = 0;
554
555 for (int d=0; d<dim; ++d)
556 {
557 for (int i=0; i<D1D[d]; ++i)
558 {
559 maxw[d] = std::max(maxw[d], maxDD[d][i] - minDD[d][i] + 1);
560 }
561 }
562
563 cdofs.SetSize(maxw[0], maxw[1], maxw[2]);
564
565 // Compute sparsity of the sparse matrix
566 smati = Memory<int>(ndof+1);
567 smati[0] = 0;
568
569 for (int dof_j=0; dof_j<ndof; ++dof_j)
570 {
571 const int jdz = dof_j / (D1D[0] * D1D[1]);
572 const int jdy = (dof_j - (jdz * D1D[0] * D1D[1])) / D1D[0];
573 const int jdx = dof_j - (jdz * D1D[0] * D1D[1]) - (jdy * D1D[0]);
574
575 MFEM_VERIFY(jdx + (D1D[0] * (jdy + (D1D[1] * jdz))) == dof_j, "");
576
577 const int jd[3] = {jdx, jdy, jdz};
578 int ndd = 1;
579 for (int i=0; i<dim; ++i)
580 {
581 nd[i] = maxDD[i][jd[i]] - minDD[i][jd[i]] + 1;
582 ndd *= nd[i];
583 }
584
585 smati[dof_j + 1] = smati[dof_j] + ndd;
586 nnz += ndd;
587 }
588
589 smatj = Memory<int>(nnz);
590 smata = Memory<real_t>(nnz);
591
592 for (int i=0; i<nnz; ++i)
593 {
594 smatj[i] = -1;
595 smata[i] = 0.0;
596 }
597
598 for (int dof_j=0; dof_j<ndof; ++dof_j)
599 {
600 const int jdz = dof_j / (D1D[0] * D1D[1]);
601 const int jdy = (dof_j - (jdz * D1D[0] * D1D[1])) / D1D[0];
602 const int jdx = dof_j - (jdz * D1D[0] * D1D[1]) - (jdy * D1D[0]);
603
604 const int jd[3] = {jdx, jdy, jdz};
605 for (int i=0; i<dim; ++i)
606 {
607 nd[i] = maxDD[i][jd[i]] - minDD[i][jd[i]] + 1;
608 }
609
610 for (int i=0; i<nd[0]; ++i)
611 for (int j=0; j<nd[1]; ++j)
612 for (int k=0; k<nd[2]; ++k)
613 {
614 cdofs(i,j,k) = minDD[0][jdx] + i + (D1D[0] *
615 (minDD[1][jdy] + j + (D1D[1] * (minDD[2][jdz] + k))));
616 }
617
618 for (int d=0; d<dim; ++d)
619 for (int qz = minD[2][jdz]; qz <= maxD[2][jdz]; ++qz)
620 {
621 for (int qy = minD[1][jdy]; qy <= maxD[1][jdy]; ++qy)
622 {
623 for (int qx = minD[0][jdx]; qx <= maxD[0][jdx]; ++qx)
624 {
625 grad[d](qx,qy,qz) = 0.0;
626 }
627 }
628 }
629
630 for (int qz = minD[2][jdz]; qz <= maxD[2][jdz]; ++qz)
631 {
632 const real_t wz = B[2](qz,jdz);
633 const real_t wDz = G[2](qz,jdz);
634
635 for (int qy = minD[1][jdy]; qy <= maxD[1][jdy]; ++qy)
636 {
637 const real_t wy = B[1](qy,jdy);
638 const real_t wDy = G[1](qy,jdy);
639
640 for (int qx = minD[0][jdx]; qx <= maxD[0][jdx]; ++qx)
641 {
642 const int q = qx + ((qy + (qz * Q1D[1])) * Q1D[0]);
643 const real_t O11 = qd(q,0);
644 const real_t O12 = qd(q,1);
645 const real_t O13 = qd(q,2);
646 const real_t O21 = symmetric ? O12 : qd(q,3);
647 const real_t O22 = symmetric ? qd(q,3) : qd(q,4);
648 const real_t O23 = symmetric ? qd(q,4) : qd(q,5);
649 const real_t O31 = symmetric ? O13 : qd(q,6);
650 const real_t O32 = symmetric ? O23 : qd(q,7);
651 const real_t O33 = symmetric ? qd(q,5) : qd(q,8);
652
653 const real_t wx = B[0](qx,jdx);
654 const real_t wDx = G[0](qx,jdx);
655
656 const real_t gradX = wDx * wy * wz;
657 const real_t gradY = wx * wDy * wz;
658 const real_t gradZ = wx * wy * wDz;
659
660 grad[0](qx,qy,qz) = (O11*gradX)+(O12*gradY)+(O13*gradZ);
661 grad[1](qx,qy,qz) = (O21*gradX)+(O22*gradY)+(O23*gradZ);
662 grad[2](qx,qy,qz) = (O31*gradX)+(O32*gradY)+(O33*gradZ);
663 }
664 }
665 }
666
667 for (int qz = minD[2][jdz]; qz <= maxD[2][jdz]; ++qz)
668 {
669 for (int dy = minDD[1][jdy]; dy <= maxDD[1][jdy]; ++dy)
670 {
671 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
672 {
673 for (int d=0; d<dim; ++d)
674 {
675 gradDXY(dx,dy,d) = 0.0;
676 }
677 }
678 }
679 for (int qy = minD[1][jdy]; qy <= maxD[1][jdy]; ++qy)
680 {
681 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
682 {
683 for (int d=0; d<dim; ++d)
684 {
685 gradDX(dx,d) = 0.0;
686 }
687 }
688 for (int qx = minD[0][jdx]; qx <= maxD[0][jdx]; ++qx)
689 {
690 const real_t gX = grad[0](qx,qy,qz);
691 const real_t gY = grad[1](qx,qy,qz);
692 const real_t gZ = grad[2](qx,qy,qz);
693 for (int dx = minQ[0][qx]; dx <= maxQ[0][qx]; ++dx)
694 {
695 const real_t wx = B[0](qx,dx);
696 const real_t wDx = G[0](qx,dx);
697 gradDX(dx,0) += gX * wDx;
698 gradDX(dx,1) += gY * wx;
699 gradDX(dx,2) += gZ * wx;
700 }
701 }
702 for (int dy = minQ[1][qy]; dy <= maxQ[1][qy]; ++dy)
703 {
704 const real_t wy = B[1](qy,dy);
705 const real_t wDy = G[1](qy,dy);
706 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
707 {
708 gradDXY(dx,dy,0) += gradDX(dx,0) * wy;
709 gradDXY(dx,dy,1) += gradDX(dx,1) * wDy;
710 gradDXY(dx,dy,2) += gradDX(dx,2) * wy;
711 }
712 }
713 }
714 for (int dz = minQ[2][qz]; dz <= maxQ[2][qz]; ++dz)
715 {
716 const real_t wz = B[2](qz,dz);
717 const real_t wDz = G[2](qz,dz);
718 for (int dy = minDD[1][jdy]; dy <= maxDD[1][jdy]; ++dy)
719 {
720 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
721 {
722 const real_t v = (gradDXY(dx,dy,0) * wz) +
723 (gradDXY(dx,dy,1) * wz) +
724 (gradDXY(dx,dy,2) * wDz);
725
726 const int loc = dx - minDD[0][jd[0]] + (nd[0] * (dy - minDD[1][jd[1]] +
727 (nd[1] * (dz - minDD[2][jd[2]]))));
728
729 const int odof = cdofs(dx - minDD[0][jd[0]],
730 dy - minDD[1][jd[1]],
731 dz - minDD[2][jd[2]]);
732
733 const int m = smati[dof_j] + loc;
734 MFEM_ASSERT(smatj[m] == odof || smatj[m] == -1, "");
735 smatj[m] = odof;
736 smata[m] += v;
737 } // dx
738 } // dy
739 } // dz
740 } // qz
741 } // dof_j
742
743 // Note that smat takes ownership of its input data.
744 smat = new SparseMatrix(smati, smatj, smata, ndof, ndof);
745}
746
747void DiffusionIntegrator::SetupPatchBasisData(Mesh *mesh, unsigned int patch)
748{
749 MFEM_VERIFY(pB.size() == patch && pG.size() == patch, "");
750 MFEM_VERIFY(pQ1D.size() == patch && pD1D.size() == patch, "");
751 MFEM_VERIFY(pminQ.size() == patch && pmaxQ.size() == patch, "");
752 MFEM_VERIFY(pminD.size() == patch && pmaxD.size() == patch, "");
753 MFEM_VERIFY(pminDD.size() == patch && pmaxDD.size() == patch, "");
754 MFEM_VERIFY(pir1d.size() == patch, "");
755
756 // Set basis functions and gradients for this patch
757 Array<const KnotVector*> pkv;
758 mesh->NURBSext->GetPatchKnotVectors(patch, pkv);
759 MFEM_VERIFY(pkv.Size() == dim, "");
760
761 Array<int> Q1D(dim);
762 Array<int> orders(dim);
763 Array<int> D1D(dim);
764 std::vector<Array2D<real_t>> B(dim);
765 std::vector<Array2D<real_t>> G(dim);
766 Array<const IntegrationRule*> ir1d(dim);
767
768 IntArrayVar2D minD(dim);
769 IntArrayVar2D maxD(dim);
770 IntArrayVar2D minQ(dim);
771 IntArrayVar2D maxQ(dim);
772
773 IntArrayVar2D minDD(dim);
774 IntArrayVar2D maxDD(dim);
775
776 for (int d=0; d<dim; ++d)
777 {
778 ir1d[d] = patchRules->GetPatchRule1D(patch, d);
779
780 Q1D[d] = ir1d[d]->GetNPoints();
781
782 orders[d] = pkv[d]->GetOrder();
783 D1D[d] = pkv[d]->GetNCP();
784
785 Vector shapeKV(orders[d]+1);
786 Vector dshapeKV(orders[d]+1);
787
788 B[d].SetSize(Q1D[d], D1D[d]);
789 G[d].SetSize(Q1D[d], D1D[d]);
790
791 minD[d].assign(D1D[d], Q1D[d]);
792 maxD[d].assign(D1D[d], 0);
793
794 minQ[d].assign(Q1D[d], D1D[d]);
795 maxQ[d].assign(Q1D[d], 0);
796
797 B[d] = 0.0;
798 G[d] = 0.0;
799
800 const Array<int>& knotSpan1D = patchRules->GetPatchRule1D_KnotSpan(patch, d);
801 MFEM_VERIFY(knotSpan1D.Size() == Q1D[d], "");
802
803 for (int i = 0; i < Q1D[d]; i++)
804 {
805 const IntegrationPoint &ip = ir1d[d]->IntPoint(i);
806 const int ijk = knotSpan1D[i];
807 const real_t kv0 = (*pkv[d])[orders[d] + ijk];
808 real_t kv1 = (*pkv[d])[0];
809 for (int j = orders[d] + ijk + 1; j < pkv[d]->Size(); ++j)
810 {
811 if ((*pkv[d])[j] > kv0)
812 {
813 kv1 = (*pkv[d])[j];
814 break;
815 }
816 }
817
818 MFEM_VERIFY(kv1 > kv0, "");
819
820 pkv[d]->CalcShape(shapeKV, ijk, (ip.x - kv0) / (kv1 - kv0));
821 pkv[d]->CalcDShape(dshapeKV, ijk, (ip.x - kv0) / (kv1 - kv0));
822
823 // Put shapeKV into array B storing shapes for all points.
824 // TODO: This should be based on NURBS3DFiniteElement::CalcShape and CalcDShape.
825 // For now, it works under the assumption that all NURBS weights are 1.
826 for (int j=0; j<orders[d]+1; ++j)
827 {
828 B[d](i,ijk + j) = shapeKV[j];
829 G[d](i,ijk + j) = dshapeKV[j];
830
831 minD[d][ijk + j] = std::min(minD[d][ijk + j], i);
832 maxD[d][ijk + j] = std::max(maxD[d][ijk + j], i);
833 }
834
835 minQ[d][i] = std::min(minQ[d][i], ijk);
836 maxQ[d][i] = std::max(maxQ[d][i], ijk + orders[d]);
837 }
838
839 // Determine which DOFs each DOF interacts with, in 1D.
840 minDD[d].resize(D1D[d]);
841 maxDD[d].resize(D1D[d]);
842 for (int i=0; i<D1D[d]; ++i)
843 {
844 const int qmin = minD[d][i];
845 minDD[d][i] = minQ[d][qmin];
846
847 const int qmax = maxD[d][i];
848 maxDD[d][i] = maxQ[d][qmax];
849 }
850 }
851
852 // Push patch data to global data structures
853 pB.push_back(B);
854 pG.push_back(G);
855
856 pQ1D.push_back(Q1D);
857 pD1D.push_back(D1D);
858
859 pminQ.push_back(minQ);
860 pmaxQ.push_back(maxQ);
861
862 pminD.push_back(minD);
863 pmaxD.push_back(maxD);
864
865 pminDD.push_back(minDD);
866 pmaxDD.push_back(maxDD);
867
868 pir1d.push_back(ir1d);
869}
870
871// This version uses reduced 1D quadrature rules.
872void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
873 const int patch, const FiniteElementSpace &fes, SparseMatrix*& smat)
874{
875 MFEM_VERIFY(patchRules, "patchRules must be defined");
876 dim = patchRules->GetDim();
877 const int spaceDim = dim; // TODO: generalize?
878
879 Mesh *mesh = fes.GetMesh();
880
881 if (VQ)
882 {
883 MFEM_VERIFY(VQ->GetVDim() == spaceDim,
884 "Unexpected dimension for VectorCoefficient");
885 }
886 if (MQ)
887 {
888 MFEM_VERIFY(MQ->GetWidth() == spaceDim,
889 "Unexpected width for MatrixCoefficient");
890 MFEM_VERIFY(MQ->GetHeight() == spaceDim,
891 "Unexpected height for MatrixCoefficient");
892 }
893
894#ifdef MFEM_THREAD_SAFE
895 DenseMatrix M(MQ ? spaceDim : 0);
896 Vector D(VQ ? VQ->GetVDim() : 0);
897#else
898 M.SetSize(MQ ? spaceDim : 0);
899 D.SetSize(VQ ? VQ->GetVDim() : 0);
900#endif
901
902 SetupPatchBasisData(mesh, patch);
903
904 MFEM_VERIFY(3 == dim, "Only 3D so far");
905
906 // Setup quadrature point data.
907 // For each point in patchRules, get the corresponding element and element
908 // reference point, in order to use element transformations. This requires
909 // data set up in NURBSPatchRule::SetPointToElement.
910 SetupPatchPA(patch, mesh, true);
911
912 const Array<int>& Q1D = pQ1D[patch];
913 const Array<int>& D1D = pD1D[patch];
914 const std::vector<Array2D<real_t>>& B = pB[patch];
915 const std::vector<Array2D<real_t>>& G = pG[patch];
916
917 const IntArrayVar2D& minD = pminD[patch];
918 const IntArrayVar2D& maxD = pmaxD[patch];
919 const IntArrayVar2D& minQ = pminQ[patch];
920 const IntArrayVar2D& maxQ = pmaxQ[patch];
921
922 const IntArrayVar2D& minDD = pminDD[patch];
923 const IntArrayVar2D& maxDD = pmaxDD[patch];
924
925 int ndof = D1D[0];
926 for (int d=1; d<dim; ++d)
927 {
928 ndof *= D1D[d];
929 }
930
931 auto rw = Reshape(reducedWeights.data(), numTypes, dim, numPatches);
932 auto rid = Reshape(reducedIDs.data(), numTypes, dim, numPatches);
933
934 const auto qd = Reshape(pa_data.Read(), Q1D[0]*Q1D[1]*Q1D[2],
935 (symmetric ? 6 : 9));
936
937 // NOTE: the following is adapted from PADiffusionApply3D.
938 std::vector<Array3D<real_t>> grad(dim);
939 for (int d=0; d<dim; ++d)
940 {
941 grad[d].SetSize(Q1D[0], Q1D[1], Q1D[2]);
942 }
943
944 // Note that this large 3D array, most of which is unused,
945 // seems inefficient, but testing showed that it is faster
946 // than using a set<int> to store 1D indices of the used points.
947 Array3D<bool> gradUsed;
948 gradUsed.SetSize(Q1D[0], Q1D[1], Q1D[2]);
949
950 Array3D<real_t> gradDXY(D1D[0], D1D[1], dim);
951 Array2D<real_t> gradDX(D1D[0], dim);
952
953 int nd[3];
954 Array3D<int> cdofs;
955
956 int *smati = nullptr;
957 int *smatj = nullptr;
958 real_t *smata = nullptr;
959 bool bugfound = false;
960 int nnz = 0;
961
962 Array<int> maxw(dim);
963 maxw = 0;
964
965 for (int d=0; d<dim; ++d)
966 {
967 for (int i=0; i<D1D[d]; ++i)
968 {
969 maxw[d] = std::max(maxw[d], maxDD[d][i] - minDD[d][i] + 1);
970 }
971 }
972
973 cdofs.SetSize(maxw[0], maxw[1], maxw[2]);
974
975 // Compute sparsity of the sparse matrix
976 smati = Memory<int>(ndof+1);
977 smati[0] = 0;
978
979 for (int dof_j=0; dof_j<ndof; ++dof_j)
980 {
981 const int jdz = dof_j / (D1D[0] * D1D[1]);
982 const int jdy = (dof_j - (jdz * D1D[0] * D1D[1])) / D1D[0];
983 const int jdx = dof_j - (jdz * D1D[0] * D1D[1]) - (jdy * D1D[0]);
984
985 MFEM_VERIFY(jdx + (D1D[0] * (jdy + (D1D[1] * jdz))) == dof_j, "");
986
987 const int jd[3] = {jdx, jdy, jdz};
988 int ndd = 1;
989 for (int i=0; i<dim; ++i)
990 {
991 nd[i] = maxDD[i][jd[i]] - minDD[i][jd[i]] + 1;
992 ndd *= nd[i];
993 }
994
995 smati[dof_j + 1] = smati[dof_j] + ndd;
996 nnz += ndd;
997 }
998
999 smatj = Memory<int>(nnz);
1000 smata = Memory<real_t>(nnz);
1001
1002 for (int i=0; i<nnz; ++i)
1003 {
1004 smatj[i] = -1;
1005 smata[i] = 0.0;
1006 }
1007
1008 for (int dof_j=0; dof_j<ndof; ++dof_j)
1009 {
1010 const int jdz = dof_j / (D1D[0] * D1D[1]);
1011 const int jdy = (dof_j - (jdz * D1D[0] * D1D[1])) / D1D[0];
1012 const int jdx = dof_j - (jdz * D1D[0] * D1D[1]) - (jdy * D1D[0]);
1013
1014 const int jd[3] = {jdx, jdy, jdz};
1015 for (int i=0; i<dim; ++i)
1016 {
1017 nd[i] = maxDD[i][jd[i]] - minDD[i][jd[i]] + 1;
1018 }
1019
1020 for (int i=0; i<nd[0]; ++i)
1021 for (int j=0; j<nd[1]; ++j)
1022 for (int k=0; k<nd[2]; ++k)
1023 {
1024 cdofs(i,j,k) = minDD[0][jdx] + i + (D1D[0] *
1025 (minDD[1][jdy] + j + (D1D[1] * (minDD[2][jdz] + k))));
1026 }
1027
1028 for (int qz = minD[2][jdz]; qz <= maxD[2][jdz]; ++qz)
1029 {
1030 for (int qy = minD[1][jdy]; qy <= maxD[1][jdy]; ++qy)
1031 {
1032 for (int qx = minD[0][jdx]; qx <= maxD[0][jdx]; ++qx)
1033 {
1034 gradUsed(qx,qy,qz) = false;
1035 }
1036 }
1037 }
1038
1039 for (int zquad = 0; zquad<2; ++zquad)
1040 {
1041 // Reduced quadrature in z
1042 const int nwz = static_cast<int>(rid(zquad,2,patch)[jdz].size());
1043 for (int irz=0; irz < nwz; ++irz)
1044 {
1045 const int qz = rid(zquad,2,patch)[jdz][irz] + minD[2][jdz];
1046 const real_t zw = rw(zquad,2,patch)[jdz][irz];
1047
1048 const real_t gwz = B[2](qz,jdz);
1049 const real_t gwDz = G[2](qz,jdz);
1050
1051 for (int dy = minDD[1][jdy]; dy <= maxDD[1][jdy]; ++dy)
1052 {
1053 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
1054 {
1055 for (int d=0; d<dim; ++d)
1056 {
1057 gradDXY(dx,dy,d) = 0.0;
1058 }
1059 }
1060 }
1061
1062 for (int yquad = 0; yquad<2; ++yquad)
1063 {
1064 // Reduced quadrature in y
1065 const int nwy = static_cast<int>(rid(yquad,1,patch)[jdy].size());
1066 for (int iry=0; iry < nwy; ++iry)
1067 {
1068 const int qy = rid(yquad,1,patch)[jdy][iry] + minD[1][jdy];
1069 const real_t yw = rw(yquad,1,patch)[jdy][iry];
1070
1071 const real_t gwy = B[1](qy,jdy);
1072 const real_t gwDy = G[1](qy,jdy);
1073
1074 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
1075 {
1076 for (int d=0; d<dim; ++d)
1077 {
1078 gradDX(dx,d) = 0.0;
1079 }
1080 }
1081
1082 // Reduced quadrature in x
1083 for (int xquad=0; xquad<2; ++xquad)
1084 {
1085 const int nwx = static_cast<int>(rid(xquad,0,patch)[jdx].size());
1086 for (int irx=0; irx < nwx; ++irx)
1087 {
1088 const int qx = rid(xquad,0,patch)[jdx][irx] + minD[0][jdx];
1089
1090 if (!gradUsed(qx,qy,qz))
1091 {
1092 const real_t gwx = B[0](qx,jdx);
1093 const real_t gwDx = G[0](qx,jdx);
1094
1095 const int q = qx + ((qy + (qz * Q1D[1])) * Q1D[0]);
1096 const real_t O11 = qd(q,0);
1097 const real_t O12 = qd(q,1);
1098 const real_t O13 = qd(q,2);
1099 const real_t O21 = symmetric ? O12 : qd(q,3);
1100 const real_t O22 = symmetric ? qd(q,3) : qd(q,4);
1101 const real_t O23 = symmetric ? qd(q,4) : qd(q,5);
1102 const real_t O31 = symmetric ? O13 : qd(q,6);
1103 const real_t O32 = symmetric ? O23 : qd(q,7);
1104 const real_t O33 = symmetric ? qd(q,5) : qd(q,8);
1105
1106 const real_t gradX = gwDx * gwy * gwz;
1107 const real_t gradY = gwx * gwDy * gwz;
1108 const real_t gradZ = gwx * gwy * gwDz;
1109
1110 grad[0](qx,qy,qz) = (O11*gradX)+(O12*gradY)+(O13*gradZ);
1111 grad[1](qx,qy,qz) = (O21*gradX)+(O22*gradY)+(O23*gradZ);
1112 grad[2](qx,qy,qz) = (O31*gradX)+(O32*gradY)+(O33*gradZ);
1113
1114 gradUsed(qx,qy,qz) = true;
1115 }
1116 }
1117 }
1118
1119 // 00 terms
1120 const int nw = static_cast<int>(rid(0,0,patch)[jdx].size());
1121 for (int irx=0; irx < nw; ++irx)
1122 {
1123 const int qx = rid(0,0,patch)[jdx][irx] + minD[0][jdx];
1124
1125 const real_t gY = grad[1](qx,qy,qz);
1126 const real_t gZ = grad[2](qx,qy,qz);
1127 const real_t xw = rw(0,0,patch)[jdx][irx];
1128 for (int dx = minQ[0][qx]; dx <= maxQ[0][qx]; ++dx)
1129 {
1130 const real_t wx = B[0](qx,dx);
1131 if (yquad == 1)
1132 {
1133 gradDX(dx,1) += gY * wx * xw;
1134 }
1135 if (zquad == 1)
1136 {
1137 gradDX(dx,2) += gZ * wx * xw;
1138 }
1139 }
1140 }
1141
1142 // 11 terms
1143 const int nw11 = static_cast<int>(rid(1,0,patch)[jdx].size());
1144
1145 for (int irx=0; irx < nw11; ++irx)
1146 {
1147 const int qx = rid(1,0,patch)[jdx][irx] + minD[0][jdx];
1148
1149 const real_t gX = grad[0](qx,qy,qz);
1150 const real_t xw = rw(1,0,patch)[jdx][irx];
1151 for (int dx = minQ[0][qx]; dx <= maxQ[0][qx]; ++dx)
1152 {
1153 const real_t wDx = G[0](qx,dx);
1154 gradDX(dx,0) += gX * wDx * xw;
1155 }
1156 }
1157
1158 for (int dy = minQ[1][qy]; dy <= maxQ[1][qy]; ++dy)
1159 {
1160 const real_t wy = B[1](qy,dy);
1161 const real_t wDy = G[1](qy,dy);
1162 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
1163 {
1164 if (yquad == 0)
1165 {
1166 if (zquad == 1)
1167 {
1168 gradDXY(dx,dy,2) += gradDX(dx,2) * wy * yw;
1169 }
1170 else
1171 {
1172 gradDXY(dx,dy,0) += gradDX(dx,0) * wy * yw;
1173 }
1174 }
1175 else if (zquad == 0)
1176 {
1177 gradDXY(dx,dy,1) += gradDX(dx,1) * wDy * yw;
1178 }
1179 }
1180 }
1181 } // qy
1182 } // y quadrature type
1183 for (int dz = minQ[2][qz]; dz <= maxQ[2][qz]; ++dz)
1184 {
1185 const real_t wz = B[2](qz,dz);
1186 const real_t wDz = G[2](qz,dz);
1187 for (int dy = minDD[1][jdy]; dy <= maxDD[1][jdy]; ++dy)
1188 {
1189 for (int dx = minDD[0][jdx]; dx <= maxDD[0][jdx]; ++dx)
1190 {
1191 real_t v = (zquad == 0) ? (gradDXY(dx,dy,0) * wz) +
1192 (gradDXY(dx,dy,1) * wz) : gradDXY(dx,dy,2) * wDz;
1193
1194 v *= zw;
1195
1196 const int loc = dx - minDD[0][jd[0]] + (nd[0] * (dy - minDD[1][jd[1]] +
1197 (nd[1] * (dz - minDD[2][jd[2]]))));
1198
1199 const int odof = cdofs(dx - minDD[0][jd[0]],
1200 dy - minDD[1][jd[1]],
1201 dz - minDD[2][jd[2]]);
1202
1203 const int m = smati[dof_j] + loc;
1204 if (!(smatj[m] == odof || smatj[m] == -1))
1205 {
1206 bugfound = true;
1207 }
1208
1209 smatj[m] = odof;
1210 smata[m] += v;
1211 } // dx
1212 } // dy
1213 } // dz
1214 } // qz
1215 } // zquad
1216 } // dof_j
1217
1218 MFEM_VERIFY(!bugfound, "");
1219
1220 for (int i=0; i<nnz; ++i)
1221 {
1222 if (smata[i] == 0.0)
1223 {
1224 // This prevents failure of SparseMatrix EliminateRowCol.
1225 // TODO: is there a better solution?
1226 smata[i] = 1.0e-16;
1227 }
1228 }
1229
1230 // Note that smat takes ownership of its input data.
1231 smat = new SparseMatrix(smati, smatj, smata, ndof, ndof);
1232}
1233
1234void DiffusionIntegrator::AssemblePatchMatrix(const int patch,
1235 const FiniteElementSpace &fes,
1236 SparseMatrix*& smat)
1237{
1238 if (integrationMode == PATCHWISE_REDUCED)
1239 {
1240 AssemblePatchMatrix_reducedQuadrature(patch, fes, smat);
1241 }
1242 else
1243 {
1244 AssemblePatchMatrix_fullQuadrature(patch, fes, smat);
1245 }
1246}
1247
1248}
mfem::Array2D
Dynamic 2D array using row-major layout.
Definition array.hpp:392
mfem::Array2D::NumCols
int NumCols() const
Definition array.hpp:408
mfem::Array2D::NumRows
int NumRows() const
Definition array.hpp:407
mfem::Array::Read
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition array.hpp:337
mfem::Coefficient::Eval
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
mfem::DenseMatrix
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
mfem::DiffusionIntegrator::AssemblePatchMatrix
void AssemblePatchMatrix(const int patch, const FiniteElementSpace &fes, SparseMatrix *&smat) override
Definition bilininteg_diffusion_patch.cpp:1234
mfem::DiffusionIntegrator::MQ
MatrixCoefficient * MQ
Definition bilininteg.hpp:2188
mfem::DiffusionIntegrator::VQ
VectorCoefficient * VQ
Definition bilininteg.hpp:2187
mfem::FiniteElementSpace
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
mfem::IntegrationPoint
Class for integration point with weight.
Definition intrules.hpp:35
mfem::IntegrationRule
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
mfem::IntegrationRule::GetNPoints
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
mfem::IntegrationRule::IntPoint
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
mfem::Integrator::patchRules
NURBSMeshRules * patchRules
Definition integrator.hpp:71
mfem::MatrixCoefficient::Eval
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
mfem::MatrixCoefficient::GetWidth
int GetWidth() const
Get the width of the matrix.
Definition coefficient.hpp:1072
mfem::MatrixCoefficient::GetHeight
int GetHeight() const
Get the height of the matrix.
Definition coefficient.hpp:1069
mfem::NNLSSolver::Mult
void Mult(const Vector &w, Vector &sol) const override
Compute the non-negative least squares solution to the underdetermined system.
Definition solvers.cpp:3728
mfem::NNLSSolver::SetOperator
void SetOperator(const Operator &op) override
The operator must be a DenseMatrix.
Definition solvers.cpp:3672
mfem::NURBSMeshRules::GetPatchRule1D
const IntegrationRule * GetPatchRule1D(const int patch, const int dimension) const
For tensor product rules defined on each patch by SetPatchRules1D(), return a pointer to the 1D rule ...
Definition intrules.hpp:314
mfem::NURBSMeshRules::GetPointElement
int GetPointElement(int patch, int i, int j, int k) const
For tensor product rules defined on each patch by SetPatchRules1D(), returns the index of the element...
Definition intrules.hpp:333
mfem::NURBSMeshRules::GetPatchRule1D_KnotSpan
const Array< int > & GetPatchRule1D_KnotSpan(const int patch, const int dimension) const
For tensor product rules defined on each patch by SetPatchRules1D(), returns an array of knot span in...
Definition intrules.hpp:343
mfem::NURBSMeshRules::GetIntegrationPointFrom1D
void GetIntegrationPointFrom1D(const int patch, int i, int j, int k, IntegrationPoint &ip)
For tensor product rules defined on each patch by SetPatchRules1D(), return the integration point wit...
Definition intrules.cpp:1968
mfem::Operator::NumCols
int NumCols() const
Get the number of columns (size of input) of the Operator. Synonym with Width().
Definition operator.hpp:75
mfem::SparseMatrix
Data type sparse matrix.
Definition sparsemat.hpp:51
mfem::VectorCoefficient::GetVDim
int GetVDim()
Returns dimension of the vector.
Definition coefficient.hpp:583
mfem::VectorCoefficient::Eval
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
mfem::Vector
Vector data type.
Definition vector.hpp:82
mfem::Vector::Read
virtual const real_t * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:494
mfem::Vector::Size
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
mfem::Vector::SetSize
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
mfem::Vector::Write
virtual real_t * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition vector.hpp:502
dim
int dim
Definition ex24.cpp:53
mfem::GetReducedRule
void GetReducedRule(const int nq, const int nd, Array2D< real_t > const &B, Array2D< real_t > const &G, std::vector< int > minQ, std::vector< int > maxQ, std::vector< int > minD, std::vector< int > maxD, std::vector< int > minDD, std::vector< int > maxDD, const IntegrationRule *ir, const bool zeroOrder, std::vector< Vector > &reducedWeights, std::vector< std::vector< int > > &reducedIDs)
Definition bilininteg_diffusion_patch.cpp:148
mfem::Reshape
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition dtensor.hpp:131
mfem::SetupPatch3D
void SetupPatch3D(const int Q1Dx, const int Q1Dy, const int Q1Dz, const int coeffDim, const bool symmetric, const Array< real_t > &w, const Vector &j, const Vector &c, Vector &d)
Definition bilininteg_diffusion_patch.cpp:25
mfem::real_t
float real_t
Definition config.hpp:43
p
real_t p(const Vector &x, real_t t)
Definition navier_mms.cpp:53
Generated on Wed Apr 9 2025 16:53:04 for MFEM by doxygen 1.11.0