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