Skip to content

Commit b0368fd

Browse files
cbritopachecocbritopacheco
committed
Update
1 parent 33ce8db commit b0368fd

File tree

7 files changed

+85
-90
lines changed

7 files changed

+85
-90
lines changed

‎examples/PETSc/MPI_Poisson.cpp‎

Lines changed: 27 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,14 @@
1-
#include "Rodin/Variational/ForwardDecls.h"
21
#include <boost/mpi/environment.hpp>
32
#include <boost/mpi/communicator.hpp>
43

54
#include <Rodin/MPI.h>
65
#include <Rodin/PETSc.h>
7-
#include <Rodin/Variational.h>
8-
9-
#include <Rodin/MPI/Geometry/Mesh.h>
10-
#include <Rodin/MPI/Geometry/Sharder.h>
116

12-
#include <Rodin/Geometry/BalancedCompactPartitioner.h>
13-
#include <string>
7+
#include <Rodin/Types.h>
8+
#include <Rodin/Solver.h>
9+
#include <Rodin/Assembly.h>
10+
#include <Rodin/Geometry.h>
11+
#include <Rodin/Variational.h>
1412

1513
namespace mpi = boost::mpi;
1614

@@ -19,6 +17,8 @@ using namespace Rodin::Math;
1917
using namespace Rodin::Solver;
2018
using namespace Rodin::Variational;
2119

20+
static constexpr Index ROOT_RANK = 0;
21+
2222
int main(int argc, char** argv)
2323
{
2424
mpi::environment env(argc, argv);
@@ -30,23 +30,37 @@ int main(int argc, char** argv)
3030
assert(ierr == PETSC_SUCCESS);
3131

3232
Rodin::MPI::Sharder sharder(mpi);
33-
if (world.rank() == 0)
33+
if (world.rank() == ROOT_RANK)
3434
{
3535
Geometry::LocalMesh mesh;
3636
mesh = mesh.UniformGrid(Geometry::Polytope::Type::Triangle, { 3, 3 });
3737
mesh.getConnectivity().compute(2, 2);
3838
mesh.getConnectivity().compute(1, 2);
3939
Geometry::BalancedCompactPartitioner partitioner(mesh);
4040
partitioner.partition(world.size());
41-
sharder.shard(partitioner).scatter(0);
41+
sharder.shard(partitioner).scatter(ROOT_RANK);
4242
}
4343

44-
auto mesh = sharder.gather(0);
44+
auto mesh = sharder.gather(ROOT_RANK);
4545
P1 vh(mesh);
4646

47-
// {
48-
// PETSc::Variational::GridFunction gf(vh);
49-
// }
47+
ScalarFunction f = 1;
48+
49+
{
50+
PETSc::Variational::TrialFunction u(vh);
51+
PETSc::Variational::TestFunction v(vh);
52+
53+
// Define problem
54+
Problem poisson(u, v);
55+
poisson = Integral(Grad(u), Grad(v))
56+
- Integral(f, v)
57+
+ DirichletBC(u, Zero());
58+
CG(poisson).solve();
59+
60+
// Save solution
61+
u.getSolution().save("Poisson.gf");
62+
mesh.save("Poisson.mesh");
63+
}
5064

5165
PetscFinalize();
5266
}

‎examples/PETSc/Seq_Poisson.cpp‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -31,7 +31,6 @@ int main(int argc, char** argv)
3131

3232
P1 vh(mesh);
3333

34-
// Leave the block so PETSc objects are cleaned up automatically !!
3534
{
3635
PETSc::Variational::TrialFunction u(vh);
3736
PETSc::Variational::TestFunction v(vh);

‎src/Rodin/Geometry.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,5 +18,6 @@
1818
#include "Geometry/Polytope.h"
1919
#include "Geometry/PolytopeTransformation.h"
2020
#include "Geometry/IsoparametricTransformation.h"
21+
#include "Geometry/BalancedCompactPartitioner.h"
2122

2223
#endif

‎src/Rodin/MPI.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#include "MPI/Context.h"
55
#include "MPI/Assembly.h"
66
#include "MPI/Variational.h"
7+
#include "MPI/Geometry.h"
78
#include "MPI/IO.h"
89

910
#endif

‎src/Rodin/MPI/Variational/P1/P1.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,11 +1,11 @@
11
#ifndef RODIN_MPI_VARIATIONAL_P1_P1_H
22
#define RODIN_MPI_VARIATIONAL_P1_P1_H
33

4+
#include <mpi.h>
45
#include <boost/serialization/optional.hpp>
5-
#include <mpi_proto.h>
66

7-
#include "Rodin/MPI/Variational/FiniteElementSpace.h"
87
#include "Rodin/MPI/Geometry/Mesh.h"
8+
#include "Rodin/MPI/Variational/FiniteElementSpace.h"
99

1010
#include "Rodin/Variational/P1/P1.h"
1111

@@ -243,7 +243,7 @@ namespace Rodin::Variational
243243
boost::optional<ElementType> send;
244244
if (localIdx)
245245
{
246-
const Index owner = shard.getOwner(d, *localIdx);
246+
const Index owner = shard.getOwner(d).at(*localIdx);
247247
if (owner == comm.rank())
248248
send = m_fes.getFiniteElement(*localIdx, globalIdx);
249249
}

‎src/Rodin/PETSc/Variational/GridFunction.h‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -98,7 +98,9 @@ namespace Rodin::Variational
9898
const auto& mesh = fes.getMesh();
9999
const auto& ctx = mesh.getContext();
100100
PetscErrorCode ierr;
101-
data.create(ctx.getCommunicator());
101+
const MPI_Comm comm = ctx.getCommunicator();
102+
ierr = VecCreate(comm, &data);
103+
assert(ierr == PETSC_SUCCESS);
102104
const size_t globalSize = fes.getSize();
103105
const size_t localSize = fes.getShard().getSize();
104106
fes.getOwnershipRange(m_begin, m_end);

‎src/Rodin/Variational/P1/QuadratureRule.h‎

Lines changed: 50 additions & 72 deletions
Original file line numberDiff line numberDiff line change
@@ -192,29 +192,35 @@ namespace Rodin::Variational
192192
::ScalarType>
193193
{
194194
public:
195-
using FESType = P1<Range, Mesh>;
195+
using FESType =
196+
P1<Range, Mesh>;
196197

197-
using LHSType = FunctionBase<LHSDerived>;
198+
using LHSType =
199+
FunctionBase<LHSDerived>;
198200

199201
using RHSType =
200202
ShapeFunctionBase<ShapeFunction<RHSDerived, FESType, TestSpace>, FESType, TestSpace>;
201203

202-
using LHSRangeType = typename FormLanguage::Traits<LHSType>::RangeType;
204+
using LHSRangeType =
205+
typename FormLanguage::Traits<LHSType>::RangeType;
203206

204-
using RHSRangeType = typename FormLanguage::Traits<RHSType>::RangeType;
207+
using RHSRangeType =
208+
typename FormLanguage::Traits<RHSType>::RangeType;
205209

206-
using IntegrandType = ShapeFunctionBase<Dot<LHSType, RHSType>>;
210+
using IntegrandType =
211+
ShapeFunctionBase<Dot<LHSType, RHSType>>;
207212

208-
using IntegrandRangeType = typename FormLanguage::Traits<IntegrandType>::RangeType;
213+
using IntegrandRangeType =
214+
typename FormLanguage::Traits<IntegrandType>::RangeType;
209215

210-
using ScalarType = typename FormLanguage::Traits<IntegrandType>::ScalarType;
216+
using ScalarType =
217+
typename FormLanguage::Traits<IntegrandType>::ScalarType;
211218

212-
using Parent = LinearFormIntegratorBase<ScalarType>;
219+
using Parent =
220+
LinearFormIntegratorBase<ScalarType>;
213221

214222
static_assert(std::is_same_v<LHSRangeType, RHSRangeType>);
215223

216-
using RangeType = Range;
217-
218224
constexpr
219225
QuadratureRule(const IntegrandType& integrand)
220226
: Parent(integrand.getLeaf()),
@@ -259,15 +265,25 @@ namespace Rodin::Variational
259265
{
260266
m_polytope = polytope;
261267
const size_t d = polytope.getDimension();
262-
const Index idx = polytope.getIndex();
263268
const auto& geometry = polytope.getGeometry();
264269
const auto& integrand = getIntegrand().getDerived();
265270
const auto& f = integrand.getLHS();
266271
const auto& fes = integrand.getFiniteElementSpace();
267-
const auto& fe = fes.getFiniteElement(d, idx);
268-
269272
const bool recompute = !m_set || m_geometry != geometry;
270-
273+
P1Element<RHSRangeType> fe(geometry);
274+
if constexpr (std::is_same_v<RHSRangeType, ScalarType>)
275+
{
276+
fe = P1Element<RHSRangeType>(geometry);
277+
}
278+
else if constexpr (std::is_same_v<RHSRangeType, Math::Vector<ScalarType>>)
279+
{
280+
// static_assert((std::is_same_v<RHSRangeType, void> || std::is_same_v<ScalarType, void>) && std::is_same_v<RHSRangeType, ScalarType>);
281+
fe = P1Element<RHSRangeType>(d, geometry);
282+
}
283+
else
284+
{
285+
assert(false);
286+
}
271287
if (recompute)
272288
{
273289
m_set = true;
@@ -282,16 +298,12 @@ namespace Rodin::Variational
282298
fe.getBasis(local)(m_basis[local], m_qf->getPoint(0));
283299
m_dot.resize(fe.getCount());
284300
}
285-
286-
static thread_local RangeType s_v;
301+
static thread_local RHSRangeType s_v;
287302
auto& p = *m_p;
288303
p.setPolytope(polytope);
289-
290304
f(s_v, p);
291-
292305
for (size_t local = 0; local < fe.getCount(); local++)
293306
m_dot[local] = Math::dot(s_v, m_basis[local]);
294-
295307
return *this;
296308
}
297309

@@ -315,7 +327,7 @@ namespace Rodin::Variational
315327
Real m_weight;
316328
Real m_distortion;
317329

318-
std::vector<RangeType> m_basis;
330+
std::vector<RHSRangeType> m_basis;
319331
std::vector<ScalarType> m_dot;
320332

321333
bool m_set;
@@ -684,8 +696,7 @@ namespace Rodin::Variational
684696
m_p(std::move(other.m_p)),
685697
m_weight(std::move(other.m_weight)),
686698
m_distortion(std::move(other.m_distortion)),
687-
m_grad1(std::move(other.m_grad1)),
688-
m_grad2(std::move(other.m_grad2)),
699+
m_grad(std::move(other.m_grad)),
689700
m_matrix(std::move(other.m_matrix))
690701
{}
691702

@@ -714,66 +725,33 @@ namespace Rodin::Variational
714725
m_weight = m_qf->getWeight(0);
715726
m_distortion = m_p->getDistortion();
716727
const size_t d = polytope.getDimension();
717-
const Index idx = polytope.getIndex();
718728
const auto& integrand = getIntegrand();
719729
const auto& lhs = integrand.getLHS();
720730
const auto& rhs = integrand.getRHS();
721731
const auto& trialfes = lhs.getFiniteElementSpace();
722732
const auto& testfes = rhs.getFiniteElementSpace();
723733
const auto& rc = m_qf->getPoint(0);
724734
const auto& p = *m_p;
725-
if (trialfes == testfes)
735+
assert(trialfes == testfes);
736+
Math::SpatialVector<ScalarType> grad(d);
737+
const P1Element<LHSRange> fe(geometry);
738+
m_matrix.resize(fe.getCount(), fe.getCount());
739+
m_grad.resize(fe.getCount());
740+
for (size_t local = 0; local < fe.getCount(); local++)
726741
{
727-
Math::SpatialVector<ScalarType> grad(d);
728-
const auto& fe = trialfes.getFiniteElement(d, idx);
729-
m_matrix.resize(fe.getCount(), fe.getCount());
730-
731-
m_grad1.resize(fe.getCount());
732-
for (size_t local = 0; local < fe.getCount(); local++)
733-
{
734-
const auto basis = fe.getBasis(local);
735-
for (size_t j = 0; j < d; j++)
736-
grad(j) = basis.template getDerivative<1>(j)(rc);
737-
m_grad1[local] = p.getJacobianInverse().transpose() * grad;
738-
}
739-
for (size_t i = 0; i < fe.getCount(); i++)
740-
m_matrix(i, i) = m_grad1[i].squaredNorm();
741-
for (size_t i = 0; i < fe.getCount(); i++)
742-
{
743-
for (size_t j = 0; j < i; j++)
744-
m_matrix(i, j) = Math::dot(m_grad1[j], m_grad1[i]);
745-
}
746-
m_matrix.template triangularView<Eigen::Upper>() = m_matrix.adjoint();
742+
const auto basis = fe.getBasis(local);
743+
for (size_t j = 0; j < d; j++)
744+
grad(j) = basis.template getDerivative<1>(j)(rc);
745+
m_grad[local] = p.getJacobianInverse().transpose() * grad;
747746
}
748-
else
747+
for (size_t i = 0; i < fe.getCount(); i++)
748+
m_matrix(i, i) = m_grad[i].squaredNorm();
749+
for (size_t i = 0; i < fe.getCount(); i++)
749750
{
750-
Math::SpatialVector<ScalarType> grad1(d), grad2(d);
751-
const auto& trialfe = lhs.getFiniteElementSpace().getFiniteElement(d, idx);
752-
m_grad1.resize(trialfe.getCount());
753-
for (size_t i = 0; i < trialfe.getCount(); i++)
754-
{
755-
const auto basis = trialfe.getBasis(i);
756-
for (size_t local = 0; local < d; local++)
757-
grad1(local) = basis.template getDerivative<1>(local)(rc);
758-
m_grad1[i] = p.getJacobianInverse().transpose() * grad1;
759-
}
760-
761-
const auto& testfe = rhs.getFiniteElementSpace().getFiniteElement(d, idx);
762-
m_grad2.resize(testfe.getCount());
763-
for (size_t i = 0; i < testfe.getCount(); i++)
764-
{
765-
const auto basis = testfe.getBasis(i);
766-
for (size_t local = 0; local < d; local++)
767-
grad2(local) = basis.template getDerivative<1>(local)(rc);
768-
m_grad2[i] = p.getJacobianInverse().transpose() * grad2;
769-
}
770-
771-
for (size_t i = 0; i < testfe.getCount(); i++)
772-
{
773-
for (size_t j = 0; j < trialfe.getCount(); j++)
774-
m_matrix(i, j) = Math::dot(m_grad1[j], m_grad2[i]);
775-
}
751+
for (size_t j = 0; j < i; j++)
752+
m_matrix(i, j) = Math::dot(m_grad[j], m_grad[i]);
776753
}
754+
m_matrix.template triangularView<Eigen::Upper>() = m_matrix.adjoint();
777755
}
778756
return *this;
779757
}
@@ -796,7 +774,7 @@ namespace Rodin::Variational
796774

797775
Real m_weight;
798776
Real m_distortion;
799-
std::vector<Math::SpatialVector<ScalarType>> m_grad1, m_grad2;
777+
std::vector<Math::SpatialVector<ScalarType>> m_grad;
800778

801779
Math::Matrix<ScalarType> m_matrix;
802780

0 commit comments

Comments
 (0)