Skip to content

Commit

Permalink
CI
Browse files Browse the repository at this point in the history
  • Loading branch information
@cbritopacheco
cbritopacheco committed Feb 11, 2025
1 parent a2cbe3a commit 819c769
Showing 1 changed file with 102 additions and 102 deletions.
204 changes: 102 additions & 102 deletions examples/PDEs/Darcy.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,127 +22,127 @@ constexpr Real pi = Math::Constants::pi();

int main(int argc, char** argv)
{
// Load mesh
Mesh mesh;
// mesh.load("Star.mesh");
mesh = mesh.UniformGrid(Polytope::Type::Triangle, { 32, 32 });
mesh.scale(1. / (31));
mesh.getConnectivity().compute(1, 2);

// Functions
P1 vh(mesh);
P1 ph(mesh, 2);

TrialFunction u(vh);
TestFunction w1(vh);

TrialFunction v(ph);
TestFunction w2(ph);

RealFunction f1 =
[](const Geometry::Point& p)
{
return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y());
};

VectorFunction f2 = VectorFunction{
[](const Geometry::Point& p)
{
return 2 * pi * (
cos(2 * pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
+ cos(2 * pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x()))
- 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y());
},
[](const Geometry::Point& p)
{
return -2 * pi * (
cos(2 * pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
+ cos(2 * pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x()))
+ 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y());
}
};
// // Load mesh
// Mesh mesh;
// // mesh.load("Star.mesh");
// mesh = mesh.UniformGrid(Polytope::Type::Triangle, { 32, 32 });
// mesh.scale(1. / (31));
// mesh.getConnectivity().compute(1, 2);

// // Functions
// P1 vh(mesh);
// P1 ph(mesh, 2);

// TrialFunction u(vh);
// TestFunction w1(vh);

// TrialFunction v(ph);
// TestFunction w2(ph);

// RealFunction f1 =
// [](const Geometry::Point& p)
// {
// return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y())
// + 2 * pi * sin(pi * p.x()) * cos(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// - 2 * pi * sin(pi * p.y()) * cos(pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x());
// return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y());
// };

// VectorFunction f2 = VectorFunction{
// [](const Geometry::Point& p)
// {
// return -2 * pi * (
// return 2 * pi * (
// cos(2 * pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// + cos(2 * pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x()))
// + 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// + pi * cos(pi * p.x()) * sin(pi * p.y());
// - 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y());
// },
// [](const Geometry::Point& p)
// {
// return 2 * pi * (
// return -2 * pi * (
// cos(2 * pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// + cos(2 * pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x()))
// + 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// + pi * sin(pi * p.x()) * cos(pi * p.y());
// + 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y());
// }
// };

// // RealFunction f1 =
// // [](const Geometry::Point& p)
// // {
// // return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y())
// // + 2 * pi * sin(pi * p.x()) * cos(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// // - 2 * pi * sin(pi * p.y()) * cos(pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x());
// // };

// // VectorFunction f2 = VectorFunction{
// // [](const Geometry::Point& p)
// // {
// // return -2 * pi * (
// // cos(2 * pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// // + cos(2 * pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x()))
// // + 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// // + pi * cos(pi * p.x()) * sin(pi * p.y());
// // },
// // [](const Geometry::Point& p)
// // {
// // return 2 * pi * (
// // cos(2 * pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// // + cos(2 * pi * p.y()) * sin(pi * p.x()) * sin(pi * p.x()))
// // + 4 * pi * pi * sin(pi * p.x()) * sin(pi * p.x()) * sin(pi * p.y()) * sin(pi * p.y())
// // + pi * sin(pi * p.x()) * cos(pi * p.y());
// // }
// // };

// auto n = BoundaryNormal(mesh);

// Problem darcy(u, v, w1, w2);
// darcy = Integral(Grad(u), Grad(w1))
// // - Integral(Div(v), w1)
// - Integral(f1, w1)
// + Integral(Jacobian(v), Jacobian(w2))
// // + Integral(u, Div(w2))
// - Integral(f2, w2)
// + DirichletBC(u, Zero())
// + DirichletBC(v, Zero(2))
// ;

// UMFPack(darcy).solve();

// mesh.save("th.mesh");
// u.getSolution().save("u.gf");
// v.getSolution().save("v.gf");

// GridFunction uEx(vh);
// uEx = [](const Geometry::Point& p) { return sin(M_PI * p.x()) * sin(M_PI * p.y()); };
// uEx.save("uEx.gf");

// GridFunction vEx(ph);
// vEx = VectorFunction{
// [](const Geometry::Point& p)
// {
// return sin(M_PI * p.x()) * sin(M_PI * p.y());
// },
// [](const Geometry::Point& p)
// {
// return sin(M_PI * p.x()) * sin(M_PI * p.y());
// }
// };
// vEx.save("vEx.gf");


// VectorFunction f3 = VectorFunction{
// [](const Geometry::Point& p)
// {
// return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y());
// },
// [](const Geometry::Point& p)
// {
// return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y());
// }
// };

auto n = BoundaryNormal(mesh);

Problem darcy(u, v, w1, w2);
darcy = Integral(Grad(u), Grad(w1))
// - Integral(Div(v), w1)
- Integral(f1, w1)
+ Integral(Jacobian(v), Jacobian(w2))
// + Integral(u, Div(w2))
- Integral(f2, w2)
+ DirichletBC(u, Zero())
+ DirichletBC(v, Zero(2))
;

UMFPack(darcy).solve();

mesh.save("th.mesh");
u.getSolution().save("u.gf");
v.getSolution().save("v.gf");

GridFunction uEx(vh);
uEx = [](const Geometry::Point& p) { return sin(M_PI * p.x()) * sin(M_PI * p.y()); };
uEx.save("uEx.gf");

GridFunction vEx(ph);
vEx = VectorFunction{
[](const Geometry::Point& p)
{
return sin(M_PI * p.x()) * sin(M_PI * p.y());
},
[](const Geometry::Point& p)
{
return sin(M_PI * p.x()) * sin(M_PI * p.y());
}
};
vEx.save("vEx.gf");


VectorFunction f3 = VectorFunction{
[](const Geometry::Point& p)
{
return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y());
},
[](const Geometry::Point& p)
{
return 2 * pi * pi * sin(pi * p.x()) * sin(pi * p.y());
}
};

TrialFunction v2(ph);
Problem darcy2(v2, w2);
darcy2 = Integral(Jacobian(v2), Jacobian(w2))
- Integral(f3, w2)
+ DirichletBC(v2, Zero(2))
;
// TrialFunction v2(ph);
// Problem darcy2(v2, w2);
// darcy2 = Integral(Jacobian(v2), Jacobian(w2))
// - Integral(f3, w2)
// + DirichletBC(v2, Zero(2))
// ;


return 0;
Expand Down

0 comments on commit 819c769

Please sign in to comment.