#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# Aufgabe 25 (b): Änderung in "1.3 Dirichlet boundary conditions" 
# (Vorlage auf der NGS-Py-Homepage)

import netgen.gui
from ngsolve import *
from netgen.geom2d import unit_square
import numpy as np

mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
Draw(mesh)
mesh.GetBoundaries()


# Dirichlet-Daten auf dem gesamten Rand von Omega
fes = H1(mesh, order=2, dirichlet="top|bottom|left|right")
#print("free dofs of fes:\n", fes.FreeDofs())

u = fes.TrialFunction()
v = fes.TestFunction()
    
# Bilinearform anpassen
a = BilinearForm(fes, symmetric=True)
a += (grad(u)*grad(v) + u*v)*dx
a.Assemble()

# Funktion g auf dem Rand
g = x*(1-x)*exp(y)
gfu = GridFunction(fes)
gfu.Set(g, BND)

    
# Linearform anpassen
f = LinearForm(fes)
f += 2*exp(y)*v*dx
f.Assemble()

r = f.vec.CreateVector()
r.data = f.vec - a.mat * gfu.vec

gfu.vec.data += a.mat.Inverse(freedofs=fes.FreeDofs()) * r  
Draw(gfu)
    
#%% exakte Lösung zum Vergleich
uex = x*(1-x)*exp(y)
# Exakte Lösung als Gitterfunktion
uexgfu = GridFunction(fes)
uexgfu.Set(uex)
Draw(uexgfu)