| 158 | |
| 159 | |
| 160 | void seissol::solver::FreeSurfaceIntegrator::initializeProjectionMatrices(unsigned maxRefinementDepth) |
| 161 | { |
| 162 | // Sub triangles |
| 163 | triRefiner.refine(maxRefinementDepth); |
| 164 | |
| 165 | const auto projectionMatrixNCols = tensor::subTriangleProjection::Shape[tensor::subTriangleProjection::index(maxRefinementDepth)][1]; |
| 166 | |
| 167 | numberOfSubTriangles = triRefiner.subTris.size(); |
| 168 | numberOfAlignedSubTriangles = tensor::subTriangleProjection::size(maxRefinementDepth) / projectionMatrixNCols; |
| 169 | |
| 170 | assert(numberOfAlignedSubTriangles * projectionMatrixNCols == tensor::subTriangleProjection::size(maxRefinementDepth)); |
| 171 | assert(numberOfSubTriangles == (1u << (2u*maxRefinementDepth))); |
| 172 | |
| 173 | const auto projectionMatrixNumberOfReals = 4 * tensor::subTriangleProjection::size(maxRefinementDepth); |
| 174 | const auto projectionMatrixMemorySize = projectionMatrixNumberOfReals * sizeof(real); |
| 175 | const auto projectionMatrixFromFaceMemoryNumberOfReals = tensor::subTriangleProjectionFromFace::size(maxRefinementDepth); |
| 176 | const auto projectionMatrixFromFaceMemorySize = projectionMatrixFromFaceMemoryNumberOfReals * sizeof(real); |
| 177 | |
| 178 | projectionMatrixMemory = |
| 179 | std::unique_ptr<real>(static_cast<real*>(seissol::memory::allocate(projectionMatrixMemorySize, Alignment))); |
| 180 | projectionMatrixFromFace = |
| 181 | std::unique_ptr<real>(static_cast<real*>(seissol::memory::allocate(projectionMatrixFromFaceMemorySize, Alignment))); |
| 182 | |
| 183 | std::fill_n(projectionMatrixMemory.get(), 0, projectionMatrixNumberOfReals); |
| 184 | std::fill_n(projectionMatrixFromFace.get(), 0, projectionMatrixFromFaceMemoryNumberOfReals); |
| 185 | |
| 186 | for (unsigned face = 0; face < 4; ++face) { |
| 187 | projectionMatrix[face] = projectionMatrixMemory.get() + face * tensor::subTriangleProjection::size(maxRefinementDepth); |
| 188 | } |
| 189 | |
| 190 | // Triangle quadrature points and weights |
| 191 | auto points = new double[numQuadraturePoints][2]; |
| 192 | auto weights = new double[numQuadraturePoints]; |
| 193 | // TODO(SW): Use the same quadrature rule, which is used for Dynamic Rupture |
| 194 | seissol::quadrature::TriangleQuadrature(points, weights, polyDegree); |
| 195 | |
| 196 | auto points3D = std::array<std::array<double, 3>, numQuadraturePoints>{}; // Points for eval of 3D basis |
| 197 | auto points2D = std::array<std::array<double, 2>, numQuadraturePoints>{}; // Points for eval of 2D basis |
| 198 | |
| 199 | // Compute projection matrices |
| 200 | for (unsigned face = 0; face < 4; ++face) { |
| 201 | for (unsigned tri = 0; tri < numberOfSubTriangles; ++tri) { |
| 202 | for (unsigned qp = 0; qp < numQuadraturePoints; ++qp) { |
| 203 | seissol::refinement::Triangle const& subTri = triRefiner.subTris[tri]; |
| 204 | const auto chiTau = std::array<double, 2>{ |
| 205 | points[qp][0] * (subTri.x[1][0] - subTri.x[0][0]) + points[qp][1] * (subTri.x[2][0] - subTri.x[0][0]) + |
| 206 | subTri.x[0][0], |
| 207 | points[qp][0] * (subTri.x[1][1] - subTri.x[0][1]) + points[qp][1] * (subTri.x[2][1] - subTri.x[0][1]) + |
| 208 | subTri.x[0][1] |
| 209 | }; |
| 210 | seissol::transformations::chiTau2XiEtaZeta(face, chiTau.data(), points3D[qp].data()); |
| 211 | points2D[qp] = chiTau; |
| 212 | } |
| 213 | computeSubTriangleAverages(projectionMatrix[face] + tri, points3D, weights); |
| 214 | if (face == 0) { |
| 215 | computeSubTriangleAveragesFromFaces(projectionMatrixFromFace.get() + tri, points2D, weights); |
| 216 | } |
| 217 | } |