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 "../linearform.hpp"
15 #include "../../general/forall.hpp"
16 #include "../../linalg/kernels.hpp"
17
18 namespace mfem
19 {
20
MFEM_REGISTER_TMOP_KERNELS(double,EnergyPA_C0_2D,const double lim_normal,const Vector & lim_dist,const Vector & c0_,const int NE,const DenseTensor & j_,const Array<double> & w_,const Array<double> & b_,const Array<double> & bld_,const Vector & x0_,const Vector & x1_,const Vector & ones,Vector & energy,const int d1d,const int q1d)21 MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
22 const double lim_normal,
23 const Vector &lim_dist,
24 const Vector &c0_,
25 const int NE,
26 const DenseTensor &j_,
27 const Array<double> &w_,
28 const Array<double> &b_,
29 const Array<double> &bld_,
30 const Vector &x0_,
31 const Vector &x1_,
32 const Vector &ones,
33 Vector &energy,
34 const int d1d,
35 const int q1d)
36 {
37 const bool const_c0 = c0_.Size() == 1;
38
39 constexpr int DIM = 2;
40 constexpr int NBZ = 1;
41
42 const int D1D = T_D1D ? T_D1D : d1d;
43 const int Q1D = T_Q1D ? T_Q1D : q1d;
44
45 const auto C0 = const_c0 ?
46 Reshape(c0_.Read(), 1, 1, 1) :
47 Reshape(c0_.Read(), Q1D, Q1D, NE);
48 const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
49 const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
50 const auto b = Reshape(b_.Read(), Q1D, D1D);
51 const auto bld = Reshape(bld_.Read(), Q1D, D1D);
52 const auto W = Reshape(w_.Read(), Q1D, Q1D);
53 const auto X0 = Reshape(x0_.Read(), D1D, D1D, DIM, NE);
54 const auto X1 = Reshape(x1_.Read(), D1D, D1D, DIM, NE);
55
56 auto E = Reshape(energy.Write(), Q1D, Q1D, NE);
57
58 MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
59 {
60 const int D1D = T_D1D ? T_D1D : d1d;
61 const int Q1D = T_Q1D ? T_Q1D : q1d;
62 constexpr int NBZ = 1;
63 constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
64 constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
65
66 MFEM_SHARED double B[MQ1*MD1];
67 MFEM_SHARED double BLD[MQ1*MD1];
68
69 MFEM_SHARED double XY[NBZ][MD1*MD1];
70 MFEM_SHARED double DQ[NBZ][MD1*MQ1];
71 MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
72
73 MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
74 MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
75 MFEM_SHARED double QQ0[2][NBZ][MQ1*MQ1];
76
77 MFEM_SHARED double XY1[2][NBZ][MD1*MD1];
78 MFEM_SHARED double DQ1[2][NBZ][MD1*MQ1];
79 MFEM_SHARED double QQ1[2][NBZ][MQ1*MQ1];
80
81 kernels::internal::LoadX<MD1,NBZ>(e,D1D,LD,XY);
82 kernels::internal::LoadX<MD1,NBZ>(e,D1D,X0,XY0);
83 kernels::internal::LoadX<MD1,NBZ>(e,D1D,X1,XY1);
84
85 kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
86 kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
87
88 kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,BLD,XY,DQ);
89 kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,BLD,DQ,QQ);
90
91 kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY0,DQ0);
92 kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ0,QQ0);
93
94 kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY1,DQ1);
95 kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ1,QQ1);
96
97 MFEM_FOREACH_THREAD(qy,y,Q1D)
98 {
99 MFEM_FOREACH_THREAD(qx,x,Q1D)
100 {
101 double ld, p0[2], p1[2];
102 const double *Jtr = &J(0,0,qx,qy,e);
103 const double detJtr = kernels::Det<2>(Jtr);
104 const double weight = W(qx,qy) * detJtr;
105 const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
106 kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ,ld);
107 kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ0,p0);
108 kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ1,p1);
109 const double dist = ld; // GetValues, default comp set to 0
110 const double id2 = 0.5 / (dist*dist);
111 const double dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
112 E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
113 }
114 }
115 });
116 return energy * ones;
117 }
118
GetLocalStateEnergyPA_C0_2D(const Vector & X) const119 double TMOP_Integrator::GetLocalStateEnergyPA_C0_2D(const Vector &X) const
120 {
121 const int N = PA.ne;
122 const int D1D = PA.maps->ndof;
123 const int Q1D = PA.maps->nqpt;
124 const int id = (D1D << 4 ) | Q1D;
125 const double ln = lim_normal;
126 const Vector &LD = PA.LD;
127 const DenseTensor &J = PA.Jtr;
128 const Array<double> &W = PA.ir->GetWeights();
129 const Array<double> &B = PA.maps->B;
130 const Array<double> &BLD = PA.maps_lim->B;
131 MFEM_VERIFY(PA.maps_lim->ndof == D1D, "");
132 MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
133 const Vector &X0 = PA.X0;
134 const Vector &C0 = PA.C0;
135 const Vector &O = PA.O;
136 Vector &E = PA.E;
137
138 MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
139 }
140
141 } // namespace mfem
142