From 9810f091ea99e7cdf8d4b28a0e6894f88908dd16 Mon Sep 17 00:00:00 2001
From: Boris Martin <boris.martin.be@gmail.com>
Date: Fri, 14 Oct 2022 11:19:57 +0200
Subject: [PATCH] working!
---
tutorials/poisson_robin_no_overlap/domain.geo | 12 ++---
tutorials/poisson_robin_no_overlap/main.cpp | 46 ++++++++-----------
2 files changed, 24 insertions(+), 34 deletions(-)
diff --git a/tutorials/poisson_robin_no_overlap/domain.geo b/tutorials/poisson_robin_no_overlap/domain.geo
index 805d938e..dbf65fdc 100644
--- a/tutorials/poisson_robin_no_overlap/domain.geo
+++ b/tutorials/poisson_robin_no_overlap/domain.geo
@@ -1,16 +1,16 @@
// Gmsh project created on Fri Oct 7 15:45:08 2022
//+
-Point(1) = {0, 0, 0, 0.05};
+Point(1) = {0, 0, 0, 0.02};
//+
-Point(2) = {1, 0, 0, 0.05};
+Point(2) = {1, 0, 0, 0.02};
//+
-Point(3) = {2, 0, 0, 0.05};
+Point(3) = {2, 0, 0, 0.02};
//+
-Point(4) = {2, 1, 0, 0.05};
+Point(4) = {2, 1, 0, 0.02};
//+
-Point(5) = {1, 1, 0, 0.05};
+Point(5) = {1, 1, 0, 0.02};
//+
-Point(6) = {0, 1, 0, 0.05};
+Point(6) = {0, 1, 0, 0.02};
//+
Line(1) = {1, 2};
//+
diff --git a/tutorials/poisson_robin_no_overlap/main.cpp b/tutorials/poisson_robin_no_overlap/main.cpp
index 0e369074..77015534 100644
--- a/tutorials/poisson_robin_no_overlap/main.cpp
+++ b/tutorials/poisson_robin_no_overlap/main.cpp
@@ -52,8 +52,8 @@ int main(int argc, char **argv)
Subdomain omega(nDom);
Interface sigma(nDom);
- omega(0) = omega0;// | interface;
- omega(1) = omega1;// | interface;
+ omega(0) = omega0;
+ omega(1) = omega1;
sigma(0, 1) = interface;
sigma(1, 0) = interface;
@@ -76,35 +76,26 @@ int main(int argc, char **argv)
InterfaceField< std::complex< double >, Form::Form0 > g("g", sigma, FunctionSpaceTypeForm0::HierarchicalH1, order);
formulation.addInterfaceField(g);
-
- // First domain
- // Laplacian(u) = -1 on the domain, and u matches g(1,0) on the interface
- formulation(0).integral(grad(dof(u(0))), grad(tf(u(0))), omega(0)|interface, gauss);
- formulation(0).integral(-1.0, tf(u(0)), omega(0)|interface, gauss);
-
- formulation(0).integral((formulation.artificialSource(-g(1, 0))), tf(u(0)), interface, gauss);
- //formulation(0).integral(dof(u(0)), tf(u(0)), interface, gauss);
+ double alpha = 1;
// First interface (What 0 sends to 1)
// Compute g(0,1) such that g(0,1) + g(1,0) = 2 u(0)
- formulation(0, 1).integral(-dof(g(0, 1)), tf(g(0, 1)), sigma(0, 1), gauss);
- formulation(0, 1).integral(-formulation.artificialSource(g(1, 0)), tf(g(0, 1)), sigma(0, 1), gauss);
- formulation(0, 1).integral((2 * u(0)), tf(g(0, 1)), sigma(0, 1), gauss);
- //formulation(0, 1).integral((2* grad(u(0))) * normal<std::complex<double>>(), tf(g(0, 1)), sigma(0, 1), gauss);
-
- // Second domain
- formulation(1).integral(grad(dof(u(1))), grad(tf(u(1))), omega(1)|interface, gauss);
- formulation(1).integral(-1.0, tf(u(1)), omega(1)|interface, gauss);
-
- formulation(1).integral(formulation.artificialSource(-g(0, 1)), tf(u(1)), interface, gauss);
- //formulation(1).integral(dof(u(1)), tf(u(1)), interface, gauss);
-
- // Second interface
- formulation(1, 0).integral(-dof(g(1, 0)), tf(g(1, 0)), sigma(1, 0), gauss);
- formulation(1, 0).integral(-formulation.artificialSource(g(0, 1)), tf(g(1, 0)), sigma(1, 0), gauss);
- formulation(1, 0).integral((2*u(1)), tf(g(1, 0)), sigma(1, 0), gauss);
- //formulation(1, 0).integral((2* grad(u(1))) * normal<std::complex<double>>(), tf(g(0, 1)), sigma(0, 1), gauss);
+ for (int i = 0; i <= 1; ++i) {
+ std::complex<double> im(0.0, 1.0);
+ auto other = 1-i;
+ // Domain
+ formulation(i).integral(grad(dof(u(i))), grad(tf(u(i))), omega(i), gauss);
+ formulation(i).integral(formulation.physicalSource(-1.0), tf(u(i)), omega(i), gauss);
+ formulation(i).integral(formulation.artificialSource(-g(other, i)), tf(u(i)), interface, gauss);
+ formulation(i).integral(alpha * (dof(u(i))), tf(u(i)), interface, gauss);
+
+ // Interface
+ formulation(i, other).integral(dof(g(i, other)), tf(g(i, other)), sigma(i, other), gauss);
+ formulation(i, other).integral(formulation.artificialSource(g(other, i)), tf(g(i, other)), sigma(i, other), gauss);
+ // It is critical that u(i) be not an artificial source (why ?)
+ formulation(i, other).integral(((-2.0 * alpha) * u(i)), tf(g(i, other)), sigma(i, other), gauss);
+ }
formulation.pre();
@@ -117,7 +108,6 @@ int main(int argc, char **argv)
if(isItMySubdomain(i)) {
save(u(i), omega(i), "u_" + std::to_string(i));
save(u(i) - trueSolution, omega(i), "err_" + std::to_string(i));
- save(g(i, 1-i), sigma(i,1-i), "g");
}
}
--
GitLab