From 98b0c1c34a2b44777ea14e0396706563f7f57e5e Mon Sep 17 00:00:00 2001
From: Jan Weidner <jw3126@gmail.com>
Date: Mon, 10 Aug 2020 11:50:45 +0200
Subject: [PATCH] add trivial.jl solving the trivial equation f = u

---
 src/trivial.jl | 39 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 39 insertions(+)
 create mode 100644 src/trivial.jl

diff --git a/src/trivial.jl b/src/trivial.jl
new file mode 100644
index 00000000..e0270337
--- /dev/null
+++ b/src/trivial.jl
@@ -0,0 +1,39 @@
+# In this tutorial we will learn
+#
+#    - How to solve the trivial equation `u = f`
+#    - How to visualize the solution in pure julia
+#
+# We want to solve `u = sin`. This equation is trivial, but it showcases how
+# the finite element machinery works. The first step is to rephrase it as
+# a variational problem:
+# ```math
+# \int u \cdot v dx = \int sin \cdot v dx
+# ```
+# for all test functions `v`.
+using Gridap
+model = CartesianDiscreteModel((0, 2π), 10) # partition the interval (0, 2π) into 10 cells
+f(pt) = sin(pt[1])
+V0 = TestFESpace(
+  reffe=:Lagrangian, order=1, valuetype=Float64,
+  conformity=:H1, model=model)
+U = TrialFESpace(V0)
+trian = Triangulation(model)
+degree = 2
+quad = CellQuadrature(trian, degree)
+A(u, v) = u ⊙ v
+b(v) = v*f
+t_Ω = AffineFETerm(A,b,trian,quad)
+op = AffineFEOperator(U,V0,t_Ω)
+u = solve(op)
+
+# Now that we have a solution, we want to know if it is any good.
+# So lets visualize it.
+using Plots
+xs = map(get_cell_coordinates(trian)) do cell
+    left, right = cell
+    1/2*(left[1] + right[1])
+end # physical coordinges of the cell centers
+q = fill([VectorValue((1/2,))], length(xs)) # reference coordinates of each cell center.
+ys = only.(evaluate(u,q)) # solution values at the cell centers
+plot(xs, ys, label="solution")
+plot!(sin, 0:0.01:2pi, label="truth")