1 // Copyright (c) 2010-2021, 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 #include "../tmop.hpp"
13 #include "tmop_pa.hpp"
14 #include "../../general/forall.hpp"
15 #include "../../linalg/kernels.hpp"
16
17 namespace mfem
18 {
19
MFEM_REGISTER_TMOP_KERNELS(void,DatcSize,const int NE,const int ncomp,const int sizeidx,const DenseMatrix & w_,const Array<double> & b_,const Vector & x_,DenseTensor & j_,const int d1d,const int q1d)20 MFEM_REGISTER_TMOP_KERNELS(void, DatcSize,
21 const int NE,
22 const int ncomp,
23 const int sizeidx,
24 const DenseMatrix &w_,
25 const Array<double> &b_,
26 const Vector &x_,
27 DenseTensor &j_,
28 const int d1d,
29 const int q1d)
30 {
31 MFEM_VERIFY(ncomp==1,"");
32 constexpr int DIM = 3;
33 const int D1D = T_D1D ? T_D1D : d1d;
34 const int Q1D = T_Q1D ? T_Q1D : q1d;
35 MFEM_VERIFY(D1D <= Q1D, "");
36
37 const auto b = Reshape(b_.Read(), Q1D, D1D);
38 const auto W = Reshape(w_.Read(), DIM,DIM);
39 const auto X = Reshape(x_.Read(), D1D, D1D, D1D, ncomp, NE);
40 auto J = Reshape(j_.Write(), DIM,DIM, Q1D,Q1D,Q1D, NE);
41
42 const double infinity = std::numeric_limits<double>::infinity();
43 MFEM_VERIFY(sizeidx == 0,"");
44 MFEM_VERIFY(MFEM_CUDA_BLOCKS==256,"");
45
46 MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
47 {
48 const int D1D = T_D1D ? T_D1D : d1d;
49 const int Q1D = T_Q1D ? T_Q1D : q1d;
50 constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
51 constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
52 constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
53
54 MFEM_SHARED double sB[MQ1*MD1];
55 MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
56 MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
57
58 ConstDeviceMatrix B(sB, D1D, Q1D);
59 DeviceCube DDD(sm0, MD1,MD1,MD1);
60 DeviceCube DDQ(sm1, MD1,MD1,MQ1);
61 DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
62 DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);
63
64 kernels::internal::LoadX(e,D1D,sizeidx,X,DDD);
65
66 double min;
67 MFEM_SHARED double min_size[MFEM_CUDA_BLOCKS];
68 DeviceTensor<3,double> M((double*)(min_size),D1D,D1D,D1D);
69 const DeviceTensor<3,const double> D((double*)(DDD+sizeidx),D1D,D1D,D1D);
70 MFEM_FOREACH_THREAD(t,x,MFEM_CUDA_BLOCKS) { min_size[t] = infinity; }
71 MFEM_SYNC_THREAD;
72 MFEM_FOREACH_THREAD(dz,z,D1D)
73 {
74 MFEM_FOREACH_THREAD(dy,y,D1D)
75 {
76 MFEM_FOREACH_THREAD(dx,x,D1D)
77 {
78 M(dx,dy,dz) = D(dx,dy,dz);
79 }
80 }
81 }
82 MFEM_SYNC_THREAD;
83 for (int wrk = MFEM_CUDA_BLOCKS >> 1; wrk > 0; wrk >>= 1)
84 {
85 MFEM_FOREACH_THREAD(t,x,MFEM_CUDA_BLOCKS)
86 { if (t < wrk) { min_size[t] = fmin(min_size[t], min_size[t+wrk]); } }
87 MFEM_SYNC_THREAD;
88 }
89 min = min_size[0];
90
91 kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,sB);
92 kernels::internal::EvalX(D1D,Q1D,B,DDD,DDQ);
93 kernels::internal::EvalY(D1D,Q1D,B,DDQ,DQQ);
94 kernels::internal::EvalZ(D1D,Q1D,B,DQQ,QQQ);
95 MFEM_FOREACH_THREAD(qx,x,Q1D)
96 {
97 MFEM_FOREACH_THREAD(qy,y,Q1D)
98 {
99 MFEM_FOREACH_THREAD(qz,z,Q1D)
100 {
101 double T;
102 kernels::internal::PullEval(qx,qy,qz,QQQ,T);
103 const double shape_par_vals = T;
104 const double size = fmax(shape_par_vals, min);
105 const double alpha = std::pow(size, 1.0/DIM);
106 for (int i = 0; i < DIM; i++)
107 {
108 for (int j = 0; j < DIM; j++)
109 {
110 J(i,j,qx,qy,qz,e) = alpha * W(i,j);
111 }
112 }
113 }
114 }
115 }
116 });
117 }
118
119 // PA.Jtr Size = (dim, dim, PA.ne*PA.nq);
ComputeAllElementTargets(const FiniteElementSpace & pa_fes,const IntegrationRule & ir,const Vector & xe,DenseTensor & Jtr) const120 void DiscreteAdaptTC::ComputeAllElementTargets(const FiniteElementSpace &pa_fes,
121 const IntegrationRule &ir,
122 const Vector &xe,
123 DenseTensor &Jtr) const
124 {
125 MFEM_VERIFY(target_type == IDEAL_SHAPE_GIVEN_SIZE ||
126 target_type == GIVEN_SHAPE_AND_SIZE,"");
127
128 MFEM_VERIFY(tspec_fesv, "No target specifications have been set.");
129 const FiniteElementSpace *fes = tspec_fesv;
130
131 // Cases that are not implemented below
132 if (skewidx != -1 ||
133 aspectratioidx != -1 ||
134 orientationidx != -1 ||
135 fes->GetMesh()->Dimension() != 3 ||
136 sizeidx == -1)
137 {
138 return ComputeAllElementTargets_Fallback(pa_fes, ir, xe, Jtr);
139 }
140
141 const Mesh *mesh = fes->GetMesh();
142 const int NE = mesh->GetNE();
143 // Quick return for empty processors:
144 if (NE == 0) { return; }
145 const int dim = mesh->Dimension();
146 MFEM_VERIFY(mesh->GetNumGeometries(dim) <= 1,
147 "mixed meshes are not supported");
148 MFEM_VERIFY(!fes->IsVariableOrder(), "variable orders are not supported");
149 const FiniteElement &fe = *fes->GetFE(0);
150 const DenseMatrix &W = Geometries.GetGeomToPerfGeomJac(fe.GetGeomType());
151 const DofToQuad::Mode mode = DofToQuad::TENSOR;
152 const DofToQuad &maps = fe.GetDofToQuad(ir, mode);
153 const Array<double> &B = maps.B;
154 const int D1D = maps.ndof;
155 const int Q1D = maps.nqpt;
156
157 Vector tspec_e;
158 const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
159 const Operator *R = fes->GetElementRestriction(ordering);
160 MFEM_VERIFY(R,"");
161 MFEM_VERIFY(R->Height() == NE*ncomp*D1D*D1D*D1D,"");
162 tspec_e.SetSize(R->Height(), Device::GetDeviceMemoryType());
163 tspec_e.UseDevice(true);
164 tspec.UseDevice(true);
165 R->Mult(tspec, tspec_e);
166 const int id = (D1D << 4 ) | Q1D;
167 MFEM_LAUNCH_TMOP_KERNEL(DatcSize,id,NE,ncomp,sizeidx,W,B,tspec_e,Jtr);
168 }
169
170 } // namespace mfem
171