From a23667f974f4b3b76f292aaf8991d15189b8b382 Mon Sep 17 00:00:00 2001
From: Nam Nguyen <nam.nguyen@rwth-aachen.de>
Date: Wed, 21 May 2025 14:34:21 +0200
Subject: [PATCH] Add FE implementations and tests for HMHF2D and RSHMHF

---
 HMHF_FE/Dim2d/__init__.py                     |   0
 HMHF_FE/Dim2d/meshing_2D.py                   |  27 ++
 HMHF_FE/Dim2d/parameters_HMHF2D.py            | 136 ++++++++
 HMHF_FE/Dim2d/solving_HMHF2D.py               | 315 ++++++++++++++++++
 HMHF_FE/Dim2d/spherical_transf.py             | 104 ++++++
 HMHF_FE/__init__.py                           |   0
 README.md                                     |  23 +-
 RSHMHF_FE/ZZ_indicator_1D.py                  |  39 +++
 RSHMHF_FE/__init__.py                         |   0
 RSHMHF_FE/adaptivity_1D.py                    | 105 ++++++
 RSHMHF_FE/marking_strategy.py                 |  27 ++
 RSHMHF_FE/meshing_1D.py                       |  67 ++++
 RSHMHF_FE/parameters_RSHMHF.py                | 133 ++++++++
 RSHMHF_FE/solving_RSHMHF.py                   | 181 ++++++++++
 Tests/__init__.py                             |   0
 Tests/conv_HMHF2D_FE.py                       | 160 +++++++++
 Tests/conv_RSHMHF_FE.py                       | 158 +++++++++
 .../input_files/HMHF_2D_PPFEM_parameters.txt  |  28 ++
 Tests/input_files/RSHMHF_parameters.txt       |  28 ++
 19 files changed, 1529 insertions(+), 2 deletions(-)
 create mode 100644 HMHF_FE/Dim2d/__init__.py
 create mode 100644 HMHF_FE/Dim2d/meshing_2D.py
 create mode 100644 HMHF_FE/Dim2d/parameters_HMHF2D.py
 create mode 100644 HMHF_FE/Dim2d/solving_HMHF2D.py
 create mode 100644 HMHF_FE/Dim2d/spherical_transf.py
 create mode 100644 HMHF_FE/__init__.py
 create mode 100644 RSHMHF_FE/ZZ_indicator_1D.py
 create mode 100644 RSHMHF_FE/__init__.py
 create mode 100644 RSHMHF_FE/adaptivity_1D.py
 create mode 100644 RSHMHF_FE/marking_strategy.py
 create mode 100644 RSHMHF_FE/meshing_1D.py
 create mode 100644 RSHMHF_FE/parameters_RSHMHF.py
 create mode 100644 RSHMHF_FE/solving_RSHMHF.py
 create mode 100644 Tests/__init__.py
 create mode 100644 Tests/conv_HMHF2D_FE.py
 create mode 100644 Tests/conv_RSHMHF_FE.py
 create mode 100644 Tests/input_files/HMHF_2D_PPFEM_parameters.txt
 create mode 100644 Tests/input_files/RSHMHF_parameters.txt

diff --git a/HMHF_FE/Dim2d/__init__.py b/HMHF_FE/Dim2d/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/HMHF_FE/Dim2d/meshing_2D.py b/HMHF_FE/Dim2d/meshing_2D.py
new file mode 100644
index 0000000..3e984be
--- /dev/null
+++ b/HMHF_FE/Dim2d/meshing_2D.py
@@ -0,0 +1,27 @@
+from netgen.geom2d import SplineGeometry
+import ngsolve as ngs
+from .parameters_HMHF2D import *
+
+def Create2DMesh(space_params):
+    domain = space_params.attributes.get("domain")
+    use_dirichlet = space_params.attributes.get("Dirichlet")
+    dirichlet_bnd_name = ""
+    if use_dirichlet:
+        dirichlet_bnd_name = space_params.attributes.get("Dirichlet_bnd_name")
+    if domain != "unit_disk":
+        raise RuntimeError(" Domain type ", domain, " is not supported! For now only unit_disk!")
+
+    geo = SplineGeometry()
+    geo.AddCircle((0,0),1,bc=dirichlet_bnd_name)
+
+    mesh_size = space_params.attributes.get("mesh_size")
+    fes_order = space_params.attributes.get("fes_order")
+
+    mesh = ngs.Mesh(geo.GenerateMesh(maxh=mesh_size))
+    #print(dirichlet_bnd_name)
+    ## use isoparametric elements for higher order
+    if fes_order > 1: 
+        #mesh.Curve(fes_order)
+        mesh.ngmesh.SecondOrder()
+    return mesh
+
diff --git a/HMHF_FE/Dim2d/parameters_HMHF2D.py b/HMHF_FE/Dim2d/parameters_HMHF2D.py
new file mode 100644
index 0000000..6e09a08
--- /dev/null
+++ b/HMHF_FE/Dim2d/parameters_HMHF2D.py
@@ -0,0 +1,136 @@
+import math
+import sys
+import fileinput
+
+class ParametersHMHF2D:
+	def __init__(self):
+		self.attributes = dict([])
+	
+	def ReadParametersFromFile(self,filename=""):
+		if len(sys.argv) < 2 and filename == "":
+			raise RuntimeError("No input file!")
+		if filename == "":
+			filename = sys.argv[1]
+		for line in fileinput.input(files=filename):
+			#print("line ", fileinput.lineno())
+			if line[0]=="#":
+				continue
+			arr_string = line.split()
+			if len(arr_string) == 0:
+				#print("Nothing in here")
+				continue
+			if len(arr_string) < 2:
+				raise RuntimeError('In line {} of the file "{}" an argument is missing!'.format(fileinput.lineno(),fileinput.filename()))
+			if len(arr_string)==2:
+				if arr_string[0] in self.attributes:
+					self.attributes[arr_string[0]]= arr_string[1]
+	def PrintParameters(self):
+		attributes_list = list(self.attributes)
+		for iterator in attributes_list:
+			print(iterator, " = ", self.attributes.get(iterator))
+
+class TimeInitialParametersHMHF2D(ParametersHMHF2D):
+	def __init__(self):
+		ParametersHMHF2D.__init__(self)
+		self.attributes.update({"tstart":0.0, "tend":1.0,\
+                "tstep":0.1,"BDF":1, "solver_type":"PPFEM"})
+
+	def ReadParametersFromFile(self,filename=""):
+		ParametersHMHF2D.ReadParametersFromFile(self,filename)
+		self.attributes["tstart"] = float(self.attributes.get("tstart"))
+		self.attributes["tend"] = float(self.attributes.get("tend"))
+		self.attributes["tstep"] = float(self.attributes.get("tstep"))
+		self.attributes["BDF"] = int(self.attributes.get("BDF"))
+		self.attributes["solver_type"] = self.attributes.get("solver_type")
+	
+	def PrintParameters(self):
+		ParametersHMHF2D.PrintParameters(self)
+
+class SpaceInitialParametersHMHF2D(ParametersHMHF2D):
+	def __init__(self):
+		ParametersHMHF2D.__init__(self)
+		self.attributes.update({"domain":"unit_disk",\
+					"mesh_size":0.1,\
+                    "Dirichlet":1, "Dirichlet_bnd_name":"circle",\
+                    "fes_order":1, "CPFEM_max_iter":100, "CPFEM_eps":1e-10})
+	def ReadParametersFromFile(self,filename=""):
+		ParametersHMHF2D.ReadParametersFromFile(self,filename)
+		self.attributes["domain"] = self.attributes.get("domain")
+		self.attributes["mesh_size"] = float(self.attributes.get("mesh_size"))
+		self.attributes["Dirichlet"] = int(self.attributes.get("Dirichlet"))
+		self.attributes["Dirichlet_bnd_name"] = self.attributes.get("Dirichlet_bnd_name")
+		self.attributes["fes_order"] = int(self.attributes.get("fes_order"))
+
+        ### parameter for fix point iteration of CPFEM
+		self.attributes["CPFEM_max_iter"] = int(self.attributes.get("CPFEM_max_iter"))
+		self.attributes["CPFEM_eps"] = float(self.attributes.get("CPFEM_eps"))
+	
+	def PrintParameters(self):
+		ParametersHMHF2D.PrintParameters(self)
+
+### ToDo: Add more output options
+class OutputParametersHMHF2D(ParametersHMHF2D):
+	def __init__(self):
+		ParametersHMHF2D.__init__(self)
+		self.attributes.update({"debug": 1, "save_mesh":"mesh.dat", \
+					"save_gfu":"gfu.dat", "save_energy":"energy.dat", \
+					"save_time":"time.dat", "save_tstep":"tstep.dat", \
+					"save_ndof":"ndof.dat", \
+					"save_temporal_error_indicator":"eta_tau.dat", \
+					"save_spacial_error_indicator":"eta_h.dat", \
+                    "save":"0","vtk":0, "vtkFilename":"", "vtkInterval":100
+					})
+
+	def ReadParametersFromFile(self,filename=""):
+		ParametersHMHF2D.ReadParametersFromFile(self,filename)
+		self.attributes["debug"] = int(self.attributes.get("debug"))
+		self.attributes["save"] = int(self.attributes.get("save"))
+		self.attributes["vtk"] = int(self.attributes.get("vtk"))
+		self.attributes["vtkFilename"] = self.attributes.get("vtkFilename")
+		self.attributes["vtkInterval"] = int(self.attributes.get("vtkInterval"))
+	
+	def PrintParameters(self):
+		ParametersHMHF2D.PrintParameters(self)
+
+class AdaptivityParametersHMHF2D(ParametersHMHF2D):
+	def __init__(self):
+		ParametersHMHF2D.__init__(self)
+		self.attributes.update({"use_adaptivity":0,"tstep_min":1e-8, \
+					"tstep_max":1e-1, "alpha_minus":0.1, \
+					"alpha_plus":2.0, "eps_err":5e-1, \
+					"eps_mesh_size":1e-12, "theta":0.5, \
+					"space_adaptive_option":"", \
+					"time_adaptive_option":"", \
+					"marking_strategy":"Doerfler", \
+					"blowup_rate":0
+					})
+
+	def ReadParametersFromFile(self,filename=""):
+		ParametersHMHF2D.ReadParametersFromFile(self,filename)
+		self.attributes["use_adaptivity"] = int(self.attributes.get("use_adaptivity"))
+		self.attributes["tstep_min"] = float(self.attributes.get("tstep_min"))
+		self.attributes["tstep_max"] = float(self.attributes.get("tstep_max"))
+		self.attributes["alpha_minus"] = float(self.attributes.get("alpha_minus"))
+		self.attributes["alpha_plus"] = float(self.attributes.get("alpha_plus"))
+		self.attributes["eps_err"] = float(self.attributes.get("eps_err"))
+		self.attributes["eps_mesh_size"] = float(self.attributes.get("eps_mesh_size"))
+		self.attributes["theta"] = float(self.attributes.get("theta"))
+	def PrintParameters(self):
+		ParametersHMHF2D.PrintParameters(self)
+
+def ReadAndPrintAllParametersFromFile(filename="",debug = 1):
+    tparams = TimeInitialParametersHMHF2D()
+    tparams.ReadParametersFromFile(filename)
+    sparams = SpaceInitialParametersHMHF2D()
+    sparams.ReadParametersFromFile(filename)
+    oparams = OutputParametersHMHF2D()
+    oparams.ReadParametersFromFile(filename)
+    aparams = AdaptivityParametersHMHF2D()
+    aparams.ReadParametersFromFile(filename)
+    if debug:
+        tparams.PrintParameters()
+        sparams.PrintParameters()
+        oparams.PrintParameters()
+        aparams.PrintParameters()
+    params = {"time":tparams,"space":sparams,"output":oparams,"adaptivity":aparams}
+    return params
diff --git a/HMHF_FE/Dim2d/solving_HMHF2D.py b/HMHF_FE/Dim2d/solving_HMHF2D.py
new file mode 100644
index 0000000..ed797eb
--- /dev/null
+++ b/HMHF_FE/Dim2d/solving_HMHF2D.py
@@ -0,0 +1,315 @@
+from ngsolve.solvers import *
+from ngsolve.comp import IntegrationRuleSpace
+
+from .parameters_HMHF2D import *
+from .meshing_2D import *
+
+import sys
+import numpy as np
+
+class SolvingHMHF2D:
+    def __init__(self,initial_cond, params):
+        self.time_init_params = params.get("time")
+        self.space_init_params = params.get("space")
+        self.output_params = params.get("output")
+        self.adaptivity_params = params.get("adaptivity")
+        self.mesh = Create2DMesh(self.space_init_params)
+        fes_order = self.space_init_params.attributes.get("fes_order")
+        dirichlet_bnd_name = self.space_init_params.attributes.get("Dirichlet_bnd_name")
+        self.solver_type_name = self.time_init_params.attributes.get("solver_type")
+        #print("Using ", self.solver_type_name, " to solve HMHF2D!")
+        if self.adaptivity_params.attributes.get("use_adaptivity"):
+            raise RuntimeError("No adaptivity implemented yet!")
+            self.fes_component = ngs.H1(self.mesh, order=fes_order,dirichlet=dirichlet_bnd_name,autoupdate=True)
+        else:
+            self.fes_component = ngs.H1(self.mesh, order=fes_order,dirichlet=dirichlet_bnd_name)
+        if self.solver_type_name == "TFEM":
+            self.fes = ngs.FESpace([self.fes_component,self.fes_component,self.fes_component,self.fes_component])
+        else:
+            #self.fes = self.fes_component**3
+            self.fes = ngs.H1(self.mesh, order=fes_order,dirichlet=dirichlet_bnd_name, dim=3)
+        #print("ndof = ", self.fes.ndof)
+
+        self.initial_cond = initial_cond
+        self.gfu = [] # store numerical solution at current timestep j [0], previous j-1 [1], j-2 [2] and so on depending on how many steps needed in bdf scheme
+
+        if self.solver_type_name != "TFEM":
+            for i in range(self.time_init_params.attributes.get("BDF")+1):
+                self.gfu.append(ngs.GridFunction(self.fes))
+                self.gfu[i].Set(self.initial_cond)
+        else:
+            for i in range(self.time_init_params.attributes.get("BDF")+1):
+                self.gfu.append(ngs.GridFunction(self.fes))
+                self.gfu[i].components[0].Set(self.initial_cond[0])
+                self.gfu[i].components[1].Set(self.initial_cond[1])
+                self.gfu[i].components[2].Set(self.initial_cond[2])
+
+            ## discrete time derivative
+            self.gfu_dot = ngs.GridFunction(self.fes)
+
+        #print("number of gfus = ",len(self.gfu), " for bdf ",self.time_init_params.attributes.get("BDF"))
+        self.gfu_extrap = ngs.GridFunction(self.fes)
+        if self.solver_type_name == "TFEM":
+            self.gfu_extrap.components[0].Set(self.initial_cond[0])
+            self.gfu_extrap.components[1].Set(self.initial_cond[1])
+            self.gfu_extrap.components[2].Set(self.initial_cond[2])
+        tstep =self.time_init_params.attributes.get("tstep") 
+        self.dt = ngs.Parameter(tstep)
+        self.bf = ngs.BilinearForm(self.fes)
+        self.lf = ngs.LinearForm(self.fes)
+
+        ### catch some cases
+        if not ( self.solver_type_name == "PPFEM" or  self.solver_type_name == "TFEM" or  self.solver_type_name == "CPFEM"):
+            raise RuntimeError("Solver type ",self.solver_type_name, " not supported! Only PPFEM, TFEM and CPFEM" )
+        if self.time_init_params.attributes.get("BDF") > 2 or self.time_init_params.attributes.get("BDF") < 1 :
+            raise RuntimeError("Only BDF 1 and 2 are supported")
+        if self.space_init_params.attributes.get("fes_order") > 2 or self.space_init_params.attributes.get("fes_order") < 1:
+            raise RuntimeError("Only P1 and P2 are supported")
+        if self.solver_type_name == "CPFEM":
+            if self.time_init_params.attributes.get("BDF") != 1:
+                raise RuntimeError("CPFEM is not supported by BDF with order higher than 1")
+            if fes_order != 1:
+                raise RuntimeError("CPFEM is not supported by finite element space with order higher than 1")
+
+        ### Initialize VTK output
+        if self.solver_type_name != "TFEM":
+            self.vtk = ngs.VTKOutput(ma = self.mesh,coefs = [self.gfu[0]],names=["gfu"],filename=self.output_params.attributes.get("vtkFilename"),subdivision = 3)
+        else:
+            self.vtk= ngs.VTKOutput(ma = self.mesh,coefs = [self.gfu[0].components[0],self.gfu[0].components[1],self.gfu[0].components[2]],names=["gfu_x","gfu_y","gfu_z"],filename=self.output_params.attributes.get("vtkFilename"),subdivision = 3)
+
+    def ConstructForms(self,use_bdf=1):
+        if self.solver_type_name == "PPFEM":
+            if use_bdf == 1:
+                self.coeff_bdf = [1,0] # [1, -1] are the actual coeff, but we homogenize the dirichlet boundary condition by solving for dgfu := gfu[0] - gfu[1]
+                self.coeff_extrap = [1]
+            elif use_bdf == 2: 
+                self.coeff_bdf =[1.5,-0.5,0.5] # [1.5, -2,0.5] same as bdf1
+                self.coeff_extrap = [2,-1]
+            else:
+                raise RuntimeError("Only BDF 1 and 2 supported!")
+            u = self.fes.TrialFunction()
+            v = self.fes.TestFunction()
+
+            ## incremental form: "unknown = gfu - gfu_prev"
+            self.bf = ngs.BilinearForm(self.fes)
+            self.bf += self.coeff_bdf[0] * u * v * ngs.dx
+            self.bf += self.dt* ngs.InnerProduct(ngs.grad(u),ngs.grad(v))* ngs.dx
+            self.bf += -self.dt * ngs.InnerProduct(ngs.grad(self.gfu_extrap),ngs.grad(self.gfu_extrap)) * u * v * ngs.dx
+            self.lf = ngs.LinearForm(self.fes)
+            self.lf += -self.dt * ngs.InnerProduct(ngs.grad(self.gfu[1]),ngs.grad(v)) * ngs.dx
+            self.lf += self.dt*ngs.InnerProduct(ngs.grad(self.gfu_extrap),ngs.grad(self.gfu_extrap)) * self.gfu[1] * v *ngs.dx
+            for i in range(1,len(self.coeff_bdf)):
+                self.lf += -self.coeff_bdf[i]*self.gfu[i]*v*ngs.dx
+
+        elif self.solver_type_name == "TFEM":
+            if use_bdf == 1:
+                self.coeff_bdf = [1,-1] 
+                self.coeff_extrap = [1]
+            elif use_bdf == 2: 
+                self.coeff_bdf =[1.5,-2,0.5]
+                self.coeff_extrap = [2,-1]
+            else:
+                raise RuntimeError("Only BDF 1 and 2 supported!")
+            u_dot_x, u_dot_y, u_dot_z,lamb = self.fes.TrialFunction()
+            v_x,v_y,v_z,w = self.fes.TestFunction()
+            self.bf = ngs.BilinearForm(self.fes)
+            self.bf +=  (u_dot_x*v_x+u_dot_y*v_y+u_dot_z*v_z)*ngs.dx
+            self.bf += 1/self.coeff_bdf[0] * self.dt * (ngs.grad(u_dot_x)*ngs.grad(v_x)+ngs.grad(u_dot_y)*ngs.grad(v_y)+ngs.grad(u_dot_z)*ngs.grad(v_z))*ngs.dx
+            self.bf += (self.gfu_extrap.components[0]*v_x+self.gfu_extrap.components[1]*v_y+self.gfu_extrap.components[2]*v_z)*lamb*ngs.dx
+            self.bf += (self.gfu_extrap.components[0]*u_dot_x+self.gfu_extrap.components[1]*u_dot_y+self.gfu_extrap.components[2]*u_dot_z)*w*ngs.dx
+
+            self.lf = ngs.LinearForm(self.fes)
+            for i in range(1,len(self.coeff_bdf)):
+                self.lf += self.coeff_bdf[i]/self.coeff_bdf[0]*(ngs.grad(self.gfu[i].components[0])*ngs.grad(v_x)+ngs.grad(self.gfu[i].components[1])*ngs.grad(v_y)+ngs.grad(self.gfu[i].components[2])*ngs.grad(v_z))*ngs.dx
+
+    def Solve(self):
+        if self.solver_type_name != "CPFEM":
+            self.Solve2(self.time_init_params.attributes.get("tstart"),self.time_init_params.attributes.get("tend"),self.time_init_params.attributes.get("BDF"))
+        else:
+            self.Solve_CrossFEM(self.time_init_params.attributes.get("tstart"),self.time_init_params.attributes.get("tend"),self.space_init_params.attributes.get("CPFEM_max_iter"),self.space_init_params.attributes.get("CPFEM_eps"))
+
+    def Solve2(self,tstart,tend,use_bdf,one_step=False):
+        t = tstart
+
+		## first steps needed to do multi-step bdf
+        bdf_counter = use_bdf
+        counter = 0
+        while bdf_counter>1 :
+            bdf_counter = bdf_counter -1
+            ## only compute one step
+            self.Solve2(t,tend,bdf_counter,True)
+            #print(" computed first time step with bdf ",bdf_counter)
+            t = t+self.dt.Get()
+            counter = counter+1
+        self.ConstructForms(use_bdf)
+        debug_mode = self.output_params.attributes.get("debug")
+        #if self.output_params.attributes.get("vtk") ==1 and use_bdf==1:
+        #    self.vtk.Do(time=0)
+        while abs(tend-t)>1e-12 and tend>t:
+            t = t+self.dt.Get()
+            if debug_mode ==1:
+                print ("\rt=", t, end="")
+                #print ("\rt=", t)
+            counter = counter + 1
+
+            ### update solution of previous steps
+            for i in range(len(self.gfu)-1,0,-1):
+                if self.solver_type_name != "TFEM":
+                    self.gfu[i].vec.data = self.gfu[i-1].vec
+                else:
+                    self.gfu[i].components[0].vec.data = self.gfu[i-1].components[0].vec
+                    self.gfu[i].components[1].vec.data = self.gfu[i-1].components[1].vec
+                    self.gfu[i].components[2].vec.data = self.gfu[i-1].components[2].vec
+
+            ### update extrapolation
+            self.gfu_extrap.vec[:]=0
+            for i in range(len(self.coeff_extrap)):
+                if self.solver_type_name != "TFEM":
+                    self.gfu_extrap.vec.data +=self.coeff_extrap[i]*self.gfu[i+1].vec
+                else:
+                    self.gfu_extrap.components[0].vec.data +=self.coeff_extrap[i]*self.gfu[i+1].components[0].vec
+                    self.gfu_extrap.components[1].vec.data +=self.coeff_extrap[i]*self.gfu[i+1].components[1].vec
+                    self.gfu_extrap.components[2].vec.data +=self.coeff_extrap[i]*self.gfu[i+1].components[2].vec
+            if self.solver_type_name == "TFEM":
+                normalized_gfu = ngs.Normalize(ngs.CoefficientFunction((self.gfu_extrap.components[0],self.gfu_extrap.components[1],self.gfu_extrap.components[2] )))
+                self.gfu_extrap.components[0].Interpolate(normalized_gfu[0])
+                self.gfu_extrap.components[1].Interpolate(normalized_gfu[1])
+                self.gfu_extrap.components[2].Interpolate(normalized_gfu[2])
+            
+            ### Solve it
+            with ngs.TaskManager():
+                self.bf.Assemble()
+                self.lf.Assemble()
+                ### TODO: put this into some postprocessing
+                if self.solver_type_name == "PPFEM":
+                    self.gfu[0].vec.data += self.bf.mat.Inverse(freedofs=self.fes.FreeDofs())*self.lf.vec
+                    normalized_gfu = ngs.Normalize(self.gfu[0])
+                    self.gfu[0].Interpolate(normalized_gfu)
+                elif self.solver_type_name == "TFEM":
+                    self.gfu_dot.vec.data = self.bf.mat.Inverse(freedofs=self.fes.FreeDofs())*self.lf.vec
+                    self.gfu[0].components[0].vec.data = 1/self.coeff_bdf[0]*self.dt.Get() * self.gfu_dot.components[0].vec
+                    self.gfu[0].components[1].vec.data = 1/self.coeff_bdf[0]*self.dt.Get() * self.gfu_dot.components[1].vec
+                    self.gfu[0].components[2].vec.data = 1/self.coeff_bdf[0]*self.dt.Get() * self.gfu_dot.components[2].vec
+                    for i in range(1,len(self.coeff_bdf)):
+                        self.gfu[0].components[0].vec.data += -self.coeff_bdf[i]/self.coeff_bdf[0]*self.gfu[i].components[0].vec
+                        self.gfu[0].components[1].vec.data += -self.coeff_bdf[i]/self.coeff_bdf[0]*self.gfu[i].components[1].vec
+                        self.gfu[0].components[2].vec.data += -self.coeff_bdf[i]/self.coeff_bdf[0]*self.gfu[i].components[2].vec
+            ### output
+            if self.output_params.attributes.get("vtk") ==1 and counter%self.output_params.attributes.get("vtkInterval") ==0:
+                self.vtk.Do(time=t) 
+
+                #print("Energy = ",ngs.Integrate(ngs.InnerProduct(ngs.grad(self.gfu[0]),ngs.grad(self.gfu[0])),self.mesh))
+            if one_step:
+                #print("only did one step")
+                break
+
+    def ChangeTimeStep(self,tstep):
+        self.dt.Set(tstep)
+
+    def Solve_CrossFEM(self,tstart,tend,max_fp_iter,epsi):
+        debug_mode = self.output_params.attributes.get("debug")
+        v = self.fes.TestFunction()
+        w = self.fes.TrialFunction()
+
+        gfu_w = ngs.GridFunction(self.fes)
+        gfu_u = ngs.GridFunction(self.fes)
+        gfu_e = ngs.GridFunction(self.fes)
+
+        gfu_w_disc_laplacian = ngs.GridFunction(self.fes)
+        gfu_e_disc_laplacian = ngs.GridFunction(self.fes)
+
+        gfu_w.Set(self.initial_cond)
+        gfu_u.Set(self.initial_cond)
+
+        ## Solution in prev step
+        gfu_w_prev = ngs.GridFunction(self.fes)
+        gfu_u_prev = ngs.GridFunction(self.fes)
+
+
+        ## For reduced L2 product
+        massLumping = ngs.IntegrationRule( points = [(0,0), (1,0), (0,1)],
+                                       weights = [1/6, 1/6, 1/6] )
+
+        ## Fixed point method formulation in incremental form:
+        a_2D = ngs.BilinearForm(self.fes)
+        a_2D+= 0.5 * self.dt * ngs.Cross(w,ngs.Cross(gfu_w_prev,gfu_w_disc_laplacian)) * v * ngs.dx(intrules={ngs.TRIG:massLumping})
+        a_2D+= w * v * ngs.dx(intrules={ngs.TRIG:massLumping})
+
+        f_2D = ngs.LinearForm(self.fes)
+        f_2D+= -0.5 * self.dt * ngs.Cross(gfu_w_prev, ngs.Cross(gfu_w_prev,gfu_w_disc_laplacian)) * v * ngs.dx(intrules={ngs.TRIG:massLumping})
+        f_2D+= -gfu_w_prev * v * ngs.dx(intrules={ngs.TRIG:massLumping})
+        f_2D+= gfu_u_prev * v * ngs.dx(intrules={ngs.TRIG:massLumping})
+
+
+        ## Forms for discrete Laplacian
+        a_discrete_laplacian = ngs.BilinearForm(self.fes)
+        a_discrete_laplacian += - w * v * ngs.dx(intrules={ngs.TRIG:massLumping})
+        a_discrete_laplacian.Assemble()
+        A_discrete_laplacian_inv =  a_discrete_laplacian.mat.Inverse(freedofs=self.fes.FreeDofs())
+
+        f_discrete_laplacian = ngs.LinearForm(self.fes)
+        f_discrete_laplacian += ngs.InnerProduct(ngs.grad(gfu_w_prev),ngs.grad(v)) * ngs.dx
+
+        f_e_discrete_laplacian = ngs.LinearForm(self.fes)
+        f_e_discrete_laplacian += ngs.InnerProduct(ngs.grad(gfu_e),ngs.grad(v)) * ngs.dx
+        t = tstart
+        count_fp_iter = 0
+        while abs(tend-t) > 1e-12:
+            t = t+self.dt.Get()
+            if debug_mode==1:
+                print ("t=", t)
+
+            gfu_u_prev.vec.data = gfu_u.vec
+            gfu_w.vec.data = gfu_u.vec
+
+            while count_fp_iter < max_fp_iter:
+                with ngs.TaskManager():
+                    gfu_w_prev.vec.data = gfu_w.vec
+                    f_discrete_laplacian.Assemble()
+
+                    gfu_w_disc_laplacian.vec.data = A_discrete_laplacian_inv * f_discrete_laplacian.vec
+
+                    a_2D.Assemble()
+                    f_2D.Assemble()
+                    inv_a = a_2D.mat.Inverse(freedofs=self.fes.FreeDofs())
+
+                    gfu_e.vec.data = inv_a * f_2D.vec
+                    gfu_w.vec.data += gfu_e.vec
+
+                    f_e_discrete_laplacian.Assemble()
+                    gfu_e_disc_laplacian.vec.data = A_discrete_laplacian_inv * f_e_discrete_laplacian.vec
+
+                    r = ngs.sqrt(ngs.Integrate(ngs.InnerProduct( ngs.Cross(gfu_w,gfu_e_disc_laplacian)+ngs.Cross(gfu_e,gfu_w_disc_laplacian),ngs.Cross(gfu_w,gfu_e_disc_laplacian)+ngs.Cross(gfu_e,gfu_w_disc_laplacian) ),self.mesh) )
+                count_fp_iter = count_fp_iter +1
+                if debug_mode==1:
+                    print("  iteration step: ", count_fp_iter, "with r=",r)
+                if r < epsi :
+                    gfu_u.vec.data = 2 * gfu_w.vec - gfu_u_prev.vec
+                    if debug_mode==1:
+                        print(" number of iterations needed: ",count_fp_iter)
+                    count_fp_iter = 0
+                    break
+            if (count_fp_iter ==  max_fp_iter ):
+                raise RuntimeError("not enough fix point iteration steps")
+
+        ### set to class variable
+        self.gfu[1].vec.data = gfu_u_prev.vec
+        self.gfu[0].vec.data = gfu_u.vec
+
+def ComputeH1Error2D(gfu_1,gfu_2,mesh_ref):
+    with ngs.TaskManager():
+            l2error = ngs.sqrt(ngs.Integrate((gfu_1-gfu_2)**2*ngs.dx,mesh_ref))
+            graderror = ngs.sqrt(ngs.Integrate(ngs.InnerProduct(ngs.grad(gfu_1)-ngs.grad(gfu_2),ngs.grad(gfu_1)-ngs.grad(gfu_2))*ngs.dx,mesh_ref))
+    h1error =ngs.sqrt(graderror**2+l2error**2) 
+    #print("\nL2 error = ",l2error)
+    #print("\nH1 error = ",h1error)
+    return l2error, h1error
+
+def ComputeH1Error2DForCF(cf_1,cf_2,grad_cf_1,grad_cf_2,mesh):
+    with ngs.TaskManager():
+        l2error = ngs.sqrt(ngs.Integrate((cf_1-cf_2)**2*ngs.dx,mesh))
+        graderror = ngs.sqrt(ngs.Integrate(ngs.InnerProduct(grad_cf_1-grad_cf_2,grad_cf_1-grad_cf_2)*ngs.dx,mesh))
+    h1error =ngs.sqrt(graderror**2+l2error**2) 
+    #print("\nL2 error = ",l2error)
+    #print("\nH1 error = ",h1error)
+    return l2error, h1error
diff --git a/HMHF_FE/Dim2d/spherical_transf.py b/HMHF_FE/Dim2d/spherical_transf.py
new file mode 100644
index 0000000..d501edb
--- /dev/null
+++ b/HMHF_FE/Dim2d/spherical_transf.py
@@ -0,0 +1,104 @@
+from netgen.geom2d import SplineGeometry
+import ngsolve as ngs
+from scipy.sparse import lil_matrix
+from scipy.sparse.linalg import spsolve
+import numpy as np
+
+## transforms RSHMHF to HMHF on unit disk in 2D
+## h_ref : mesh size for the mesh_2D of gfu_2D (comparable to mesh size of RSHMHF reference solution gfu_1D
+## initial_cond: 3-vectorial 2-dim initial condition of HMHF2D
+## FirstOrder = fes first order = nodal basis
+def SphericalTransform_FirstOrder(gfu_1D,mesh_1D,h_ref,initial_cond):
+    geo = SplineGeometry()
+    geo.AddCircle((0,0),1,bc="unit_circle")
+    mesh_2D = ngs.Mesh(geo.GenerateMesh(maxh=h_ref))
+    fes_2D = ngs.H1(mesh_2D,order=1,dirichlet=geo.GetBCName(1))
+    gfu_2D_1 = ngs.GridFunction(fes_2D)
+    gfu_2D_1.Set(initial_cond[0])
+    gfu_2D_2 = ngs.GridFunction(fes_2D)
+    gfu_2D_2.Set(initial_cond[1])
+    gfu_2D_3 = ngs.GridFunction(fes_2D)
+    gfu_2D_3.Set(initial_cond[2])
+    j=0
+    for p in mesh_2D.ngmesh.Points():
+        r_p = ngs.sqrt(p[0]**2 + p[1]**2);
+        if abs(r_p -1) < 1e-12: # dont change anything at the boundary
+          j = j+1
+          continue
+        gfu_2D_1.vec[j] = p[0]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0)))
+        gfu_2D_2.vec[j] = p[1]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0)))
+        gfu_2D_3.vec[j] = ngs.cos(gfu_1D(mesh_1D(r_p,0.0)))
+        j = j+1
+    gfu_2D = ngs.CoefficientFunction((gfu_2D_1,gfu_2D_2,gfu_2D_3))
+    
+    ## Jacobian of gfu_2D, compareable to ngs.grad(gfu_2D)
+    gradgfu_2D = ngs.CoefficientFunction((ngs.grad(gfu_2D_1)[0],ngs.grad(gfu_2D_2)[0],ngs.grad(gfu_2D_3)[0],ngs.grad(gfu_2D_1)[1],ngs.grad(gfu_2D_2)[1],ngs.grad(gfu_2D_3)[1]), dims=(2,3))
+    ## Gradient of gfu_2D, compareable to ngs.Grad(gfu_2D)
+    #gradgfu_2D = ngs.CoefficientFunction((ngs.grad(gfu_2D_1)[0],ngs.grad(gfu_2D_1)[1],ngs.grad(gfu_2D_2)[0],ngs.grad(gfu_2D_2)[1],ngs.grad(gfu_2D_3)[0],ngs.grad(gfu_2D_3)[1]), dims=(3,2))
+    return mesh_2D, fes_2D, gfu_2D, gradgfu_2D
+
+## transforms RSHMHF (short "1D") to HMHF on unit disk in 2D
+## h_ref : mesh size for the mesh_2D of gfu_2D (comparable to mesh size of RSHMHF reference solution gfu_1D
+## initial_cond: to set for the boundary
+## fes_order: should be same order as fes of gfu_1D
+def SphericalTransform(gfu_1D,mesh_1D,h_ref,initial_cond, fes_order):
+    if fes_order == 1:
+        mesh_2D, fes_2D, gfu_2D, gradgfu_2D = SphericalTransform_FirstOrder(gfu_1D,mesh_1D,h_ref,initial_cond)
+        return mesh_2D, fes_2D, gfu_2D, gradgfu_2D
+    
+    geo = SplineGeometry()
+    geo.AddCircle((0,0),1,bc="unit_circle")
+
+    mesh_2D = ngs.Mesh(geo.GenerateMesh(maxh=h_ref))
+    #mesh_2D.Curve(2)
+    mesh_2D.ngmesh.SecondOrder()
+    #fes_2D = ngs.H1(mesh_2D,order=fes_order,dirichlet=geo.GetBCName(1))
+    fes_2D = ngs.H1(mesh_2D,order=fes_order)
+    gfu_2D_1 = ngs.GridFunction(fes_2D)
+    gfu_2D_1.Set(initial_cond[0])
+    gfu_2D_2 = ngs.GridFunction(fes_2D)
+    gfu_2D_2.Set(initial_cond[1])
+    gfu_2D_3 = ngs.GridFunction(fes_2D)
+    gfu_2D_3.Set(initial_cond[2])
+    gfu_fes = ngs.GridFunction(fes_2D)
+
+    fes_ndof = fes_2D.ndof
+
+
+    A_1 = lil_matrix((fes_ndof,fes_ndof),dtype = np.float32) 
+    f_1 = np.zeros([fes_ndof,1])
+    f_2 = np.zeros([fes_ndof,1])
+    f_3 = np.zeros([fes_ndof,1])
+    x_ng_pnt = np.ndarray(shape=fes_2D.ndof)
+    y_ng_pnt = np.ndarray(shape=fes_2D.ndof)
+    j=0
+    for p in mesh_2D.ngmesh.Points():
+        x_ng_pnt[j] = p[0]
+        y_ng_pnt[j] = p[1]
+        r_p =ngs.sqrt(p[0]**2+p[1]**2)
+        f_1[j,0] = p[0]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0,0.0)))
+        f_2[j,0] = p[1]/r_p * ngs.sin(gfu_1D(mesh_1D(r_p,0.0,0.0)))
+        f_3[j,0] = ngs.cos(gfu_1D(mesh_1D(r_p,0.0,0.0)))
+        j+=1
+    pnt = mesh_2D(x_ng_pnt,y_ng_pnt)
+    for i in range(fes_ndof):
+        gfu_fes.vec[:] = 0
+        gfu_fes.vec[i] = 1
+        A_1[i] = lil_matrix(gfu_fes(pnt)).transpose()
+    A = A_1.transpose()
+    eps = 1e-12
+    for i in range(A.shape[0]):
+        row_data = A.data[i]
+        row_cols = A.rows[i]
+        A.data[i] = [val for val in row_data if abs(val) >= eps]
+        A.rows[i] = [col for val, col in zip(row_data, row_cols) if abs(val) >= eps]
+    A = A.tocsr()
+    gfu_2D_1.vec.data= spsolve(A,f_1)
+    gfu_2D_2.vec.data= spsolve(A,f_2)
+    gfu_2D_3.vec.data= spsolve(A,f_3)
+    #print("solving system successfully")
+
+    gfu_2D = ngs.CoefficientFunction((gfu_2D_1,gfu_2D_2,gfu_2D_3))
+    gradgfu_2D = ngs.CoefficientFunction((ngs.grad(gfu_2D_1)[0],ngs.grad(gfu_2D_2)[0],ngs.grad(gfu_2D_3)[0],ngs.grad(gfu_2D_1)[1],ngs.grad(gfu_2D_2)[1],ngs.grad(gfu_2D_3)[1]), dims=(2,3))
+
+    return mesh_2D, fes_2D, gfu_2D, gradgfu_2D
diff --git a/HMHF_FE/__init__.py b/HMHF_FE/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/README.md b/README.md
index eb068d7..69366de 100644
--- a/README.md
+++ b/README.md
@@ -2,9 +2,28 @@
 
 ## Installation
 
-[matlab](https://www.mathworks.com/products/matlab.html) are needed to run the code.
+[matlab](https://www.mathworks.com/products/matlab.html) (at least R2022b) and [ngsolve](https://www.ngsolve.org) (latest version) are needed to run the code.
 
 ## Usage
 
-- In `RSHMHF_FD`: Finite difference discretizations with [matlab](https://www.mathworks.com/products/matlab.html) for radially symmetric HMHF (RSHMHF) used in [arXiv:2503.03525](https://arxiv.org/abs/2503.03525)
+- `RSHMHF_FD`: BDF1/BDF2 + Finite difference discretizations with [matlab](https://www.mathworks.com/products/matlab.html) for radially symmetric HMHF (RSHMHF)
+  - "Discretization error analysis for a radially symmetric harmonic map heat flow problem" - Nguyen, Reusken (2025) [arXiv](https://arxiv.org/abs/2503.03525)
+- `RSHMHF_FE`: BDF1/BDF2 + P1/P2 Finite element discretizations with [ngsolve](https://www.ngsolve.org) for RSHMHF
+  - `convex_lin`: "Error analysis for a finite element discretization of a radially symmetric harmonic map heat flow problem" - Nguyen, Reusken (2025)
+  - `max_lin`: "A comparative study of finite element methods for a class of harmonic map heat flow problems" - Nguyen, Reusken (2025)
+- `HMHF_FE`: BDF1/BDF2 + P1/P2 Finite element discretizations wtih [ngsolve](https://www.ngsolve.org) for HMHF
+  - `Dim2d`:
+    - `PPFEM`: p.139 of "Computational Micromagnetism" - Andreas Prohl (2001) [DOI](https://doi.org/10.1007/978-3-663-09498-2)
+    - `TFEM`: "Higher-order linearly implicit full discretization of the Landau-Lifshitz-Gilbert equation" - Akrivis, Feischl, Kovács, Lubich (2021) [DOI](https://dx.doi.org/10.1090/mcom/3597)
+    - `CPFEM`: "Constraint preserving implicit finite element discretization of harmonic map flow into spheres" - Bartels, Prohl (2007) [DOI](https://dx.doi.org/10.1090/S0025-5718-07-02026-1)
+
+## Tests
+Run the following in command line to test mesh size and time step convergence for P1 and BDF1 of RSHMHF FE
+```console
+foo@bar:~$ python3 -m Tests.conv_RSHMHF_FE
+
+Run the following in command line to test mesh size and time step convergence for P1 and BDF1 of PPFEM and TFEM of HMHF FE on unit disk in dimension 2
+```console
+foo@bar:~$ python3 -m Tests.conv_HMHF2D_FE
+```
 
diff --git a/RSHMHF_FE/ZZ_indicator_1D.py b/RSHMHF_FE/ZZ_indicator_1D.py
new file mode 100644
index 0000000..db20920
--- /dev/null
+++ b/RSHMHF_FE/ZZ_indicator_1D.py
@@ -0,0 +1,39 @@
+import ngsolve as ngs
+
+# p0_fes := order 0 L2-space of mesh
+# gfu is current solution and lives in solution space H1 order 1 on mesh
+# return: piecewise gradient of gfu on mesh in p0_fes_space
+def ComputeP0Gradient(mesh,p0_fes,gfu):
+	assert p0_fes.ndof+1 == len(gfu.vec)
+	p0_grad = ngs.GridFunction(p0_fes)
+	mesh_size_element_wise = ngs.Integrate(ngs.CoefficientFunction(1)*ngs.dx,mesh,element_wise=True)
+	for i in range(0, len(p0_grad.vec)):
+		p0_grad.vec[i] = (gfu.vec[i+1] - gfu.vec[i])/mesh_size_element_wise[i]
+	return p0_grad
+
+# p1_fes := order 1 H1 space on mesh
+# p0_grad : piecewise constant gradient in order 0 L2 on mesh
+# return : ZZ post processed gradient in p1_fes on mesh
+def ComputeZZGradient(mesh,p1_fes,p0_grad):
+	assert p1_fes.ndof == len(p0_grad.vec)+1
+	zz_grad_gfu = ngs.GridFunction(p1_fes)
+	
+	#print("\n",p1_fes.ndof)
+	#print(len(zz_grad_gfu.vec))
+	for i in range(0,p1_fes.ndof):
+		if i==0:
+			zz_grad_gfu.vec[0] = p0_grad.vec[0]
+			continue
+		if i==p1_fes.ndof-1:
+			zz_grad_gfu.vec[p1_fes.ndof-1] = p0_grad.vec[p1_fes.ndof-2]
+			continue
+		zz_grad_gfu.vec[i] = 0.5*(p0_grad.vec[i] + p0_grad.vec[i-1])
+	return zz_grad_gfu
+
+def ComputeZZIndicator(mesh,p1_fes,gfu):
+	W = ngs.L2(mesh,order = 0)
+	p0_grad_gfu = ComputeP0Gradient(mesh,W,gfu)
+	zz_grad_gfu = ComputeZZGradient(mesh,p1_fes,p0_grad_gfu)
+	err_elwise = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh,element_wise = True)
+	err_total = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh)
+	return err_elwise, err_total
diff --git a/RSHMHF_FE/__init__.py b/RSHMHF_FE/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/RSHMHF_FE/adaptivity_1D.py b/RSHMHF_FE/adaptivity_1D.py
new file mode 100644
index 0000000..4a1a064
--- /dev/null
+++ b/RSHMHF_FE/adaptivity_1D.py
@@ -0,0 +1,105 @@
+from .parameters_RSHMHF import *
+from .meshing_1D import *
+
+from .ZZ_indicator_1D import *
+from .marking_strategy import *
+
+# gfu_list := [gfu at current time j, gfu at j-1, gfu at j-2,... ] in p1_fes on mesh
+# p1_fes := order 1 H1 space on mesh
+# tstep := current time step
+# tcurrent := current time
+# tend := end time
+# number_elements := number of elements in mesh
+# time_adaptivity_type: "TimeDerivativeGradient"
+# space_adaptive_type: "ZZRecovery", "Energy"
+# marking_strategy: "Doerfler", "equidistributed"
+def SpaceTimeAdaptivity(mesh,gfu_list,gfu_extrap,p1_fes,tstep,tcurrent,tend,number_elements,adapt_params,counter_tmin,stop_space_adapt = False,stop_time_adapt=False,debug = 0):
+    use_adaptivity = adapt_params.attributes.get("use_adaptivity")
+    tstep_min = adapt_params.attributes.get("tstep_min")
+    tstep_max = adapt_params.attributes.get("tstep_max")
+    alpha_plus = adapt_params.attributes.get("alpha_plus")
+    alpha_minus = adapt_params.attributes.get("alpha_minus")
+    eps_err = adapt_params.attributes.get("eps_err")
+    eps_mesh_size = adapt_params.attributes.get("eps_mesh_size")
+    theta = adapt_params.attributes.get("theta")
+    space_adaptive_option = adapt_params.attributes.get("space_adaptive_option")
+    time_adaptive_option = adapt_params.attributes.get("time_adaptive_option")
+    marking_strategy = adapt_params.attributes.get("marking_strategy")
+    new_tstep = tstep
+    new_tcurrent = tcurrent
+    new_mesh = mesh
+    new_gfu_list = gfu_list
+    new_number_elements = number_elements
+    new_gfu_extrap = gfu_extrap
+    dirichlet_bnd = mesh.GetBoundaries()[0]+ "|" + mesh.GetBoundaries()[1]
+    #print(dirichlet_bnd)
+    new_p1_fes = ngs.H1(mesh,order=1,dirichlet = dirichlet_bnd)
+    refined_something = False
+    if not use_adaptivity:
+        return new_mesh,new_p1_fes,new_gfu_list, new_gfu_extrap, new_number_elements, new_tstep, new_tcurrent,refined_something
+    time_indicator = 0.0
+    #time_indicator_elwise = []
+    space_indicator = 0.0
+    space_indicator_elwise = []
+    W = ngs.L2(mesh,order = 0)
+    if time_adaptive_option == "TimeDerivativeGradient" and not stop_time_adapt:
+        p0_grad_gfu = ComputeP0Gradient(mesh,W,gfu_list[0])
+        p0_grad_gfu_prev = ComputeP0Gradient(mesh,W,gfu_list[1])
+        time_indicator = ngs.Integrate((p0_grad_gfu-p0_grad_gfu_prev)**2 * ngs.dx, mesh)
+        #time_indicator_elwise = ngs.Integrate((p0_grad_gfu-p0_grad_gfu_prev)**2 * ngs.dx, mesh, element_wise=True)
+    elif time_adaptive_option == "TimeDerivativeGradient" and stop_time_adapt:
+        time_indicator = 0.0
+    elif time_adaptive_option != "TimeDerivativeGradient":
+        raise RuntimeError("Time indicator of type ",time_adaptive_option, " not supported")
+
+    if space_adaptive_option == "ZZRecovery" and not stop_space_adapt:
+        p0_grad_gfu = ComputeP0Gradient(mesh,W,gfu_list[0])
+        zz_grad_gfu = ComputeZZGradient(mesh,p1_fes,p0_grad_gfu)
+        space_indicator_elwise = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh,element_wise = True)
+        space_indicator = ngs.Integrate((zz_grad_gfu-p0_grad_gfu)**2*ngs.dx,mesh)
+        #space_indicator_elwise,space_indicator = ComputeZZIndicator(mesh,p1_fes,gfu_list[0])
+    elif space_adaptive_option == "ZZRecovery" and stop_space_adapt:
+        space_indicator = 0
+    elif space_adaptive_option != "ZZRecovery":
+        raise RuntimeError("Space indicator of type ",space_adaptive_option, " not supported")
+    #print("\n eps_err = ",eps_err)
+
+    if space_indicator + time_indicator < eps_err:
+        refined_something = False
+        # incrase time step if possible
+        if not stop_time_adapt:
+            new_tstep = min([tstep_max,alpha_plus*tstep,tend-tcurrent])
+            counter_tmin = 0
+    else:
+		# reject current solution
+        gfu_list[0].vec.data = gfu_list[1].vec
+        refined_something = True
+        new_tcurrent = tcurrent-tstep
+        ref_flags = [False]*number_elements
+        # if minimal time step is used too many times without decreasing error, then refine instead
+        if counter_tmin > 3:
+            time_indicator = 0
+            counter_tmin = 0
+		# check which error dominates
+        if time_indicator > space_indicator:
+            # reduce time step if possible
+            new_tstep = max([tstep_min,alpha_minus*tstep])
+            if abs(new_tstep -tstep_min)<1e-12:
+                counter_tmin = counter_tmin +1
+        else:
+            index_to_be_refined = SortErrListAndReturnElementIndex(number_elements,space_indicator_elwise,theta,space_indicator)
+            for index in index_to_be_refined:
+                ref_flags[index] = True
+            new_mesh, new_number_elements = Refine1DMesh(mesh,ref_flags)
+            new_p1_fes = ngs.H1(new_mesh,order=1,dirichlet=dirichlet_bnd)
+            new_gfu_list = []
+            for i in range(0,len(gfu_list)):
+                new_gfu_list.append(ngs.GridFunction(new_p1_fes))
+                new_gfu_list[i].Interpolate(gfu_list[i])
+            new_gfu_extrap = ngs.GridFunction(new_p1_fes)
+            new_gfu_extrap.Interpolate(gfu_extrap)
+            counter_tmin = 0
+    if debug:
+        print("\n time indicator = ",time_indicator)
+        print("\n space indicator = ",space_indicator)
+    return new_mesh,new_p1_fes,new_gfu_list, new_gfu_extrap, new_number_elements, new_tstep, new_tcurrent,refined_something, counter_tmin
diff --git a/RSHMHF_FE/marking_strategy.py b/RSHMHF_FE/marking_strategy.py
new file mode 100644
index 0000000..4f65cf2
--- /dev/null
+++ b/RSHMHF_FE/marking_strategy.py
@@ -0,0 +1,27 @@
+import ngsolve as ngs
+import numpy as np
+
+# Doerfler strategy
+def SortErrListAndReturnElementIndex(num_elem,err_elem,theta,err_total):
+    assert len(err_elem) == num_elem
+    elerr_index = [i for i in range(0,num_elem)]
+    elerr_sorted = list( sorted(err_elem,reverse=True))
+    elerr_index_sorted =list( np.argsort(err_elem)[::-1][:num_elem])
+    index_refine =[]
+    elerr_refine = []
+    #print("total error",err_total)
+    #print("err list")
+    #for i,j in zip(elerr,elerr_index):
+    #   print ("%.12f"%i,j)
+    #print("")
+    #print("sorted err list")
+    #for i,j in zip(elerr_sorted,elerr_index_sorted):
+    #   print ("%.12f"%i,j)
+    #print("")
+    while sum(elerr_refine) < theta * err_total and len(elerr_index_sorted) > 0:
+       elerr_refine.append(elerr_sorted.pop(0))
+       index_refine.append(elerr_index_sorted.pop(0))
+       #for i,j in zip(elerr_refine,index_refine):
+       #   print ("%.12f"%i,j)
+       #print("sum of elerr_refine so far ", "%.12f"%sum(elerr_refine), "; is it less than ", theta*err_total)
+    return index_refine
diff --git a/RSHMHF_FE/meshing_1D.py b/RSHMHF_FE/meshing_1D.py
new file mode 100644
index 0000000..dc0851e
--- /dev/null
+++ b/RSHMHF_FE/meshing_1D.py
@@ -0,0 +1,67 @@
+import netgen.meshing as ng
+import ngsolve as ngs
+
+def Print1DMesh(mesh):
+	current_discr_points =[p[0] for p in mesh.ngmesh.Points()]
+	print(current_discr_points)
+
+def Create1DMesh(number_elements,interval_start,interval_end,region_name_left = "left",region_name_right="right",region_name_inner = "material"):
+	m = ng.Mesh(dim=1)
+	a = interval_start
+	b = interval_end
+	interval_length = b-a
+	assert(interval_length > 0)
+	assert(number_elements > 0)
+	pnums = []
+	for i in range(0, number_elements+1):
+		pnums.append (m.Add (ng.MeshPoint (ng.Pnt(interval_start + interval_length * i/number_elements, 0, 0))))
+
+	idx = m.AddRegion(region_name_inner, dim=1)
+	for i in range(0,number_elements):
+		m.Add (ng.Element1D ([pnums[i],pnums[i+1]], index=idx))
+
+	idx_left = m.AddRegion(region_name_left, dim=0)
+	idx_right = m.AddRegion(region_name_right, dim=0)
+
+	m.Add (ng.Element0D (pnums[0], index=idx_left))
+	m.Add (ng.Element0D (pnums[number_elements], index=idx_right))
+
+	mesh = ngs.Mesh(m)
+
+	return mesh
+
+def Create1DMeshFromPoints(discr_points,region_name_left = "left",region_name_right="right",region_name_inner = "material"):
+	pnums = []
+	m = ng.Mesh(dim=1)
+	for i in range(0, len(discr_points)):
+		pnums.append (m.Add (ng.MeshPoint (ng.Pnt(discr_points[i], 0, 0))))
+	idx = m.AddRegion(region_name_inner, dim=1)
+	for i in range(0,len(discr_points)-1):
+		m.Add (ng.Element1D ([pnums[i],pnums[i+1]], index=idx))
+
+	idx_left = m.AddRegion(region_name_left, dim=0)
+	idx_right = m.AddRegion(region_name_right, dim=0)
+	m.Add (ng.Element0D (pnums[0], index=idx_left))
+	m.Add (ng.Element0D (pnums[len(discr_points)-1], index=idx_right))
+	mesh = ngs.Mesh(m)
+	return mesh
+
+# returns refined mesh and number of elements of the refined mesh
+# number_elements + 1 = number_discretization in 1D
+def Refine1DMesh(mesh,refinement_flags):
+	current_discr_points =[p[0] for p in mesh.ngmesh.Points()]
+	number_discretization_pts = len(current_discr_points)
+	assert number_discretization_pts == len(current_discr_points)
+	assert number_discretization_pts == len(refinement_flags)+1
+	assert len(mesh.GetMaterials()) == 1
+	assert len(mesh.GetBoundaries()) == 2
+	ref_discr_points = current_discr_points
+	nmb_added_new_points = 0
+	for i in range(0,len(refinement_flags)):
+		if refinement_flags[i]:
+			ref_discr_points.insert(i+1+nmb_added_new_points, 0.5*(current_discr_points[i+nmb_added_new_points]+current_discr_points[i+1+nmb_added_new_points]))
+			nmb_added_new_points = nmb_added_new_points +1
+	#print(ref_discr_points)
+	number_elements_refined_mesh = len(ref_discr_points)-1
+	mesh_new = Create1DMeshFromPoints(ref_discr_points, mesh.GetBoundaries()[0],mesh.GetBoundaries()[1],mesh.GetMaterials()[0])
+	return mesh_new, number_discretization_pts + nmb_added_new_points-1
diff --git a/RSHMHF_FE/parameters_RSHMHF.py b/RSHMHF_FE/parameters_RSHMHF.py
new file mode 100644
index 0000000..37ccc54
--- /dev/null
+++ b/RSHMHF_FE/parameters_RSHMHF.py
@@ -0,0 +1,133 @@
+import math
+import sys
+import fileinput
+
+class ParametersRSHMHF:
+	def __init__(self):
+		self.attributes = dict([])
+	
+	def ReadParametersFromFile(self,filename=""):
+		if len(sys.argv) < 2 and filename == "":
+			raise RuntimeError("No input file!")
+		if filename == "":
+			filename = sys.argv[1]
+		for line in fileinput.input(files=filename):
+			#print("line ", fileinput.lineno())
+			if line[0]=="#":
+				continue
+			arr_string = line.split()
+			if len(arr_string) == 0:
+				#print("Nothing in here")
+				continue
+			if len(arr_string) < 2:
+				raise RuntimeError('In line {} of the file "{}" an argument is missing!'.format(fileinput.lineno(),fileinput.filename()))
+			if len(arr_string)==2:
+				if arr_string[0] in self.attributes:
+					self.attributes[arr_string[0]]= arr_string[1]
+	def PrintParameters(self):
+		attributes_list = list(self.attributes)
+		for iterator in attributes_list:
+			print(iterator, " = ", self.attributes.get(iterator))
+
+class TimeInitialParametersRSHMHF(ParametersRSHMHF):
+	def __init__(self):
+		ParametersRSHMHF.__init__(self)
+		self.attributes.update({"tstart":0.0, "tend":1.0,\
+                "tstep":0.1,"BDF":1,"solver_type":"convex_lin"})
+
+	def ReadParametersFromFile(self,filename=""):
+		ParametersRSHMHF.ReadParametersFromFile(self,filename)
+		self.attributes["tstart"] = float(self.attributes.get("tstart"))
+		self.attributes["tend"] = float(self.attributes.get("tend"))
+		self.attributes["tstep"] = float(self.attributes.get("tstep"))
+		self.attributes["BDF"] = int(self.attributes.get("BDF"))
+		self.attributes["solver_type"] = self.attributes.get("solver_type")
+	
+	def PrintParameters(self):
+		ParametersRSHMHF.PrintParameters(self)
+
+class SpaceInitialParametersRSHMHF(ParametersRSHMHF):
+	def __init__(self):
+		ParametersRSHMHF.__init__(self)
+		self.attributes.update({"interval_start":0.0, "interval_end":1.0, \
+					"number_elements":10,"refinement_level":0,\
+					"refinement_until_this_point":1.0, \
+                    "Dirichlet":1, "fes_order":1})
+	def ReadParametersFromFile(self,filename=""):
+		ParametersRSHMHF.ReadParametersFromFile(self,filename)
+		self.attributes["interval_start"] = float(self.attributes.get("interval_start"))
+		self.attributes["interval_end"] = float(self.attributes.get("interval_end"))
+		self.attributes["number_elements"] = int(self.attributes.get("number_elements"))
+		self.attributes["refinement_level"] = int(self.attributes.get("refinement_level"))
+		self.attributes["refinement_until_this_point"] \
+			 = float(self.attributes.get("refinement_until_this_point"))
+		self.attributes["Dirichlet"] = int(self.attributes.get("Dirichlet"))
+		self.attributes["fes_order"] = int(self.attributes.get("fes_order"))
+	
+	def PrintParameters(self):
+		ParametersRSHMHF.PrintParameters(self)
+
+### ToDo: Add output options
+class OutputParametersRSHMHF(ParametersRSHMHF):
+	def __init__(self):
+		ParametersRSHMHF.__init__(self)
+		self.attributes.update({"debug": 1, "save_mesh":"mesh.dat", \
+					"save_gfu":"gfu.dat", "save_energy":"energy.dat", \
+					"save_time":"time.dat", "save_tstep":"tstep.dat", \
+					"save_ndof":"ndof.dat", \
+					"save_temporal_error_indicator":"eta_tau.dat", \
+					"save_spacial_error_indicator":"eta_h.dat", \
+					"save":"0"
+					})
+
+	def ReadParametersFromFile(self,filename=""):
+		ParametersRSHMHF.ReadParametersFromFile(self,filename)
+		self.attributes["debug"] = int(self.attributes.get("debug"))
+		self.attributes["save"] = int(self.attributes.get("save"))
+	
+	def PrintParameters(self):
+		ParametersRSHMHF.PrintParameters(self)
+
+class AdaptivityParametersRSHMHF(ParametersRSHMHF):
+	def __init__(self):
+		ParametersRSHMHF.__init__(self)
+		self.attributes.update({"use_adaptivity":0,"tstep_min":1e-8, \
+					"tstep_max":1e-1, "alpha_minus":0.1, \
+					"alpha_plus":2.0, "eps_err":5e-1, \
+					"eps_mesh_size":1e-12, "theta":0.5, \
+					"space_adaptive_option":"", \
+					"time_adaptive_option":"", \
+					"marking_strategy":"Doerfler", \
+					"blowup_rate":0
+					})
+
+	def ReadParametersFromFile(self,filename=""):
+		ParametersRSHMHF.ReadParametersFromFile(self,filename)
+		self.attributes["use_adaptivity"] = int(self.attributes.get("use_adaptivity"))
+		self.attributes["tstep_min"] = float(self.attributes.get("tstep_min"))
+		self.attributes["tstep_max"] = float(self.attributes.get("tstep_max"))
+		self.attributes["alpha_minus"] = float(self.attributes.get("alpha_minus"))
+		self.attributes["alpha_plus"] = float(self.attributes.get("alpha_plus"))
+		self.attributes["eps_err"] = float(self.attributes.get("eps_err"))
+		self.attributes["eps_mesh_size"] = float(self.attributes.get("eps_mesh_size"))
+		self.attributes["theta"] = float(self.attributes.get("theta"))
+		self.attributes["blowup_rate"] = float(self.attributes.get("blowup_rate"))
+	def PrintParameters(self):
+		ParametersRSHMHF.PrintParameters(self)
+
+def ReadAndPrintAllParametersFromFile(filename="",debug = 1):
+    tparams = TimeInitialParametersRSHMHF()
+    tparams.ReadParametersFromFile(filename)
+    sparams = SpaceInitialParametersRSHMHF()
+    sparams.ReadParametersFromFile(filename)
+    oparams = OutputParametersRSHMHF()
+    oparams.ReadParametersFromFile(filename)
+    aparams = AdaptivityParametersRSHMHF()
+    aparams.ReadParametersFromFile(filename)
+    if debug:
+        tparams.PrintParameters()
+        sparams.PrintParameters()
+        oparams.PrintParameters()
+        aparams.PrintParameters()
+    return tparams,sparams,oparams,aparams
+
diff --git a/RSHMHF_FE/solving_RSHMHF.py b/RSHMHF_FE/solving_RSHMHF.py
new file mode 100644
index 0000000..f6199a8
--- /dev/null
+++ b/RSHMHF_FE/solving_RSHMHF.py
@@ -0,0 +1,181 @@
+from ngsolve.solvers import *
+
+from .adaptivity_1D import *
+
+import sys
+import numpy as np
+
+class Blowuprate():
+    def __init__(self):
+        self.dgfudx = []
+        self.time = []
+        self.value = 0.0
+
+
+class SolvingRSHMHF:
+    def __init__(self, mesh, fes, initial_cond, params):
+        self.mesh = mesh
+        self.fes = fes
+        self.initial_cond = initial_cond
+        self.time_init_params = params.get("time")
+        ### ToDo: add some output options
+        self.output_params = params.get("output")
+        self.adaptivity_params = params.get("adaptivity")
+        self.gfu = [] # store numerical solution at current timestep j [0], previous j-1 [1], j-2 [2] and so on depending on how many steps needed in bdf scheme
+        for i in range(self.time_init_params.attributes.get("BDF")+1):
+        	self.gfu.append(ngs.GridFunction(self.fes))
+        	self.gfu[i].Set(self.initial_cond)
+        	#self.gfu[i].vec[0] = 0.0
+        #print("number of gfus = ",len(self.gfu), " for bdf ",self.time_init_params.attributes.get("BDF"))
+        self.gfu_extrap = ngs.GridFunction(self.fes)
+        tstep =self.time_init_params.attributes.get("tstep") 
+        self.dt = ngs.Parameter(tstep)
+        self.bf = ngs.BilinearForm(self.fes)
+        self.lf = ngs.LinearForm(self.fes)
+        self.blowuprate = Blowuprate()
+
+    def ConstructForms(self,use_bdf=1):
+        if use_bdf == 1:
+            coeff_bdf = [1,0] # [1, -1] are the actual coeff, but we homogenize the dirichlet boundary condition by solving for dgfu := gfu[0] - gfu[1]
+            coeff_extrap = [1]
+        elif use_bdf == 2: 
+            coeff_bdf =[1.5,-0.5,0.5] # [1.5, -2,0.5] same as bdf1
+            coeff_extrap = [2,-1]
+        else:
+            raise RuntimeError("Only BDF 1 and 2 supported!")
+        self.gfu_extrap.vec[:]=0
+        for i in range(len(coeff_extrap)):
+            self.gfu_extrap.vec.data +=coeff_extrap[i]*self.gfu[i+1].vec
+        u = self.fes.TrialFunction()
+        v = self.fes.TestFunction()
+        self.bf = ngs.BilinearForm(self.fes)
+        self.lf = ngs.LinearForm(self.fes)
+        if self.time_init_params.attributes.get("solver_type")=="convex_lin":
+            self.bf += coeff_bdf[0] * u * v *ngs.CoefficientFunction((ngs.x)) *ngs.dx
+            self.bf += self.dt* ngs.grad(u) * ngs.grad(v)*ngs.CoefficientFunction((ngs.x))* ngs.dx
+            self.bf += self.dt* u * v*ngs.CoefficientFunction((1/(ngs.x))) *ngs.dx
+
+            self.lf += - self.dt * ngs.grad(self.gfu[1])* ngs.grad(v)*ngs.CoefficientFunction((ngs.x)) * ngs.dx
+            self.lf += - self.dt* self.gfu[1] * v*ngs.CoefficientFunction((1/(ngs.x))) *ngs.dx
+            self.lf += self.dt*  ngs.CoefficientFunction((1/(ngs.x))) *self.gfu_extrap* v *ngs.dx
+            self.lf += - self.dt * ngs.CoefficientFunction((1/(ngs.x))) * ngs.sin(2*self.gfu_extrap)/(2)  * v * ngs.dx
+
+            for i in range(1,len(coeff_bdf)):
+                self.lf += -coeff_bdf[i]*self.gfu[i]*v*ngs.CoefficientFunction((ngs.x))*ngs.dx
+        elif self.time_init_params.attributes.get("solver_type")=="max_lin":
+            self.bf += coeff_bdf[0] * u * v * ngs.dx
+            self.bf += self.dt* ngs.grad(u) * ngs.grad(v)* ngs.dx
+            self.bf += -self.dt * ngs.CoefficientFunction((1/ngs.x)) * ngs.grad(u) * v* ngs.dx
+            self.bf += self.dt * ngs.CoefficientFunction((1/(ngs.x*ngs.x))) * ngs.sin(2*self.gfu_extrap)/(2*self.gfu_extrap)*u * v* ngs.dx
+            self.lf += - self.dt * ngs.grad(self.gfu[1])* ngs.grad(v) * ngs.dx
+            self.lf += self.dt * ngs.CoefficientFunction((1/ngs.x)) * ngs.grad(self.gfu[1]) * v *ngs.dx
+            self.lf += -self.dt * ngs.CoefficientFunction((1/(ngs.x*ngs.x))) * ngs.sin(2*self.gfu_extrap)/(2*self.gfu_extrap) * self.gfu[1] * v * ngs.dx
+            for i in range(1,len(coeff_bdf)):
+                self.lf += -coeff_bdf[i]*self.gfu[i]*v*ngs.dx
+        else:
+            raise RuntimeError("solver type ",self.time_init_params.attributes.get("solver_type"), " not supported. Use convex_lin or max_lin!")
+
+    def Solve(self):
+        self.Solve2(self.time_init_params.attributes.get("tstart"),self.time_init_params.attributes.get("tend"),self.time_init_params.attributes.get("BDF"))
+
+    def Solve2(self,tstart,tend,use_bdf,one_step=False):
+        if use_bdf == 1:
+            coeff_bdf = [1.0,0.0] # [1, -1] are the actual coeff, but we homogenize the dirichlet boundary condition by solving for dgfu := gfu[0] - gfu[1]
+            coeff_extrap = [1]
+        elif use_bdf == 2: 
+            coeff_bdf =[1.5,-0.5,0.5] # [1.5, -2,0.5] same as bdf1
+            coeff_extrap = [2,-1]
+        else:
+            raise RuntimeError("Only BDF 1 and 2 supported!")
+        new_tstep = self.dt.Get()
+        number_elements = self.fes.ndof -1 # linear fes
+        new_mesh = self.mesh
+        new_fes = self.fes
+        new_gfu = self.gfu
+        new_gfu_extrap = self.gfu_extrap
+        t = tstart
+
+		## first steps needed to do multi-step bdf
+        bdf_counter = use_bdf
+        while bdf_counter>1 :
+            bdf_counter = bdf_counter -1
+            self.Solve2(t,tend,bdf_counter,True)
+            t = t+self.dt.Get()
+        self.ConstructForms(use_bdf)
+        debug_mode = self.output_params.attributes.get("debug")
+        use_adaptivity = self.adaptivity_params.attributes.get("use_adaptivity")
+        counter = 0
+        tsteps_arr = []
+        counter_tmin = 0
+        while abs(tend-t)>1e-12 and tend>t:
+            t = t+self.dt.Get()
+            counter = counter + 1
+            if debug_mode:
+                print("\rt=",t,end="")
+            
+            for i in range(len(self.gfu)-1,0,-1):
+            	self.gfu[i].vec.data = self.gfu[i-1].vec
+            
+            self.gfu_extrap.vec[:]=0
+            for i in range(len(coeff_extrap)):
+            	self.gfu_extrap.vec.data +=coeff_extrap[i]*self.gfu[i+1].vec
+            
+            with ngs.TaskManager():
+            	self.bf.Assemble()
+            	self.lf.Assemble()
+            	self.gfu[0].vec.data += self.bf.mat.Inverse(freedofs=self.fes.FreeDofs())*self.lf.vec
+
+            if use_adaptivity == 1:
+                self.blowuprate.value = self.adaptivity_params.attributes.get("blowup_rate")
+                if use_bdf > 1:
+                    raise RuntimeError("Only BDF1 supported in adaptive mode!")
+                new_mesh,new_fes,new_gfu,new_gfu_extrap,number_elements,new_tstep,new_tcurrent,refined_something,counter_tmin = SpaceTimeAdaptivity(self.mesh,self.gfu,self.gfu_extrap,self.fes,self.dt.Get(),t,tend,len(self.gfu[0].vec)-1,self.adaptivity_params,counter_tmin)
+                mesh_size_el = ngs.Integrate(ngs.CoefficientFunction(1)*ngs.dx,self.mesh,element_wise=True)
+                min_mesh_size = min(mesh_size_el)
+                if min_mesh_size < self.adaptivity_params.attributes.get("eps_mesh_size"):
+                    print("\nReach minimal mesh size of ", min_mesh_size)
+                    break
+                self.dt.Set(new_tstep)
+                if refined_something:
+                    self.mesh = new_mesh
+                    self.fes = ngs.H1(self.mesh,order=1,dirichlet="left|right")
+                    self.gfu = []
+                    for i in range(0,len(new_gfu)):
+                        self.gfu.append(ngs.GridFunction(self.fes))
+                        self.gfu[i].Interpolate(new_gfu[i])
+                    self.gfu_extrap = ngs.GridFunction(self.fes)
+                    self.gfu_extrap.Interpolate(new_gfu_extrap)
+                    t = new_tcurrent
+                    self.ConstructForms(use_bdf)
+                    counter = counter - 1
+                    continue
+            if self.blowuprate.value>0:
+                self.blowuprate.dgfudx.append((self.gfu[0])(self.mesh(self.mesh.ngmesh.Points()[2][0],0.0,0.0))/self.mesh.ngmesh.Points()[2][0] )
+                self.blowuprate.time.append(t)
+            if debug_mode and use_adaptivity:
+            	print("\n ndof = ", self.fes.ndof)
+            	print("\n timestep = ",self.dt.Get())
+            if one_step:
+            	break
+        if debug_mode and use_adaptivity:
+            print(" Number of time steps needed = ",counter)
+            print("\n ndof = ", self.fes.ndof)
+
+def ComputeH1Error(gfu_1,gfu_2,mesh_ref):
+    with ngs.TaskManager():
+        l2error = ngs.sqrt(ngs.Integrate((gfu_1-gfu_2)**2*ngs.dx,mesh_ref))
+        graderror = ngs.sqrt(ngs.Integrate((ngs.grad(gfu_1)-ngs.grad(gfu_2))**2*ngs.dx,mesh_ref))
+    h1error =ngs.sqrt(graderror**2+l2error**2) 
+    #print("L2 error = ",l2error)
+    #print("H1 error = ",h1error)
+    return l2error, h1error
+
+def ComputeH1RAndEnergyError(gfu_1,gfu_2,mesh_ref):
+    with ngs.TaskManager():
+        l2rerror = ngs.sqrt(ngs.Integrate((gfu_1-gfu_2)**2*ngs.CoefficientFunction((ngs.x))*ngs.dx,mesh_ref))
+        h1rerror = ngs.sqrt(ngs.Integrate(((ngs.grad(gfu_1)-ngs.grad(gfu_2))**2*ngs.CoefficientFunction((ngs.x))+(gfu_1-gfu_2)**2*ngs.CoefficientFunction((1/ngs.x)))*ngs.dx,mesh_ref))
+        energyerror = ngs.sqrt(ngs.Integrate(((ngs.grad(gfu_1)-ngs.grad(gfu_2))**2*ngs.CoefficientFunction((ngs.x))+ngs.sin((gfu_1-gfu_2))**2*ngs.CoefficientFunction((1/ngs.x)))*ngs.dx,mesh_ref))
+    #print("L2R error = ",l2rerror)
+    #print("H1R error = ",h1rerror)
+    #print("Energy error = ",energyerror)
+    return l2rerror, h1rerror, energyerror
diff --git a/Tests/__init__.py b/Tests/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/Tests/conv_HMHF2D_FE.py b/Tests/conv_HMHF2D_FE.py
new file mode 100644
index 0000000..1573a0a
--- /dev/null
+++ b/Tests/conv_HMHF2D_FE.py
@@ -0,0 +1,160 @@
+## Add project root to sys.path
+import sys
+import os
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+
+from RSHMHF_FE import parameters_RSHMHF
+from RSHMHF_FE.solving_RSHMHF import *
+from HMHF_FE.Dim2d  import parameters_HMHF2D
+from HMHF_FE.Dim2d.solving_HMHF2D  import *
+from HMHF_FE.Dim2d.spherical_transf import *
+import os
+
+print("############### Compute reference solution from RSHMHF ####################")
+
+## initial condition for RSHMHF
+## requirement: u_0(x=0) = 0
+u0 = ngs.CoefficientFunction((0.5*math.pi*ngs.x*ngs.x))
+
+## parameter input file
+this_dir = os.path.dirname(__file__)
+input_file = os.path.join(this_dir, "input_files/RSHMHF_parameters.txt")
+time_init_params_ref,space_init_params_ref,output_params_ref, adaptivity_params_ref =parameters_RSHMHF.ReadAndPrintAllParametersFromFile(input_file,0)
+time_init_params_ref.attributes.get("solver_type") 
+number_elements_ref = space_init_params_ref.attributes.get("number_elements")
+use_bdf_ref = time_init_params_ref.attributes.get("BDF")
+dt_ref = time_init_params_ref.attributes.get("tstep")
+fes_order = space_init_params_ref.attributes.get("fes_order")
+
+print("use bdf ", use_bdf_ref)
+print("time steps = ",dt_ref)
+print("ref number elements = ",number_elements_ref)
+print("fes_order = ",fes_order)
+print("use adaptivity = ",adaptivity_params_ref.attributes.get("use_adaptivity"))
+a = space_init_params_ref.attributes.get("interval_start")
+b = space_init_params_ref.attributes.get("interval_end")
+mesh_ref = Create1DMesh(number_elements_ref,a,b)
+
+use_dirichlet = space_init_params_ref.attributes.get("Dirichlet")
+if use_dirichlet:
+    ## by default: "left" and "right"
+	dirichlet_str = mesh_ref.GetBoundaries()[0] + "|" + mesh_ref.GetBoundaries()[1]
+	#print(dirichlet_str)
+	V_ref = ngs.H1(mesh_ref,order=fes_order,dirichlet=dirichlet_str)
+	#print(V_ref.FreeDofs())
+else:
+	raise RuntimeError("Only with Dirichlet boundary condition!")
+
+print("ndof = ",V_ref.ndof)
+
+params_ref = {"time":time_init_params_ref,"output":output_params_ref,"adaptivity":adaptivity_params_ref}
+params_ref["time"].attributes["solver_type"] = "max_lin"
+solver_RSHMHF_ref = SolvingRSHMHF(mesh_ref,V_ref,u0,params_ref)
+solver_RSHMHF_ref.Solve()
+
+print("\n############### Compute solution HMHF 2D ####################")
+
+## input file
+this_dir = os.path.dirname(__file__)
+input_file = os.path.join(this_dir, "input_files/HMHF_2D_PPFEM_parameters.txt")
+params = parameters_HMHF2D.ReadAndPrintAllParametersFromFile(input_file,0)
+
+time_init_params_ref = params.get("time")
+space_init_params_ref = params.get("space")
+output_params_ref = params.get("output")
+adaptivity_params_ref = params.get("adaptivity")
+fes_order = space_init_params_ref.attributes.get("fes_order")
+
+print("fes order = ",fes_order)
+print("domain = ", space_init_params_ref.attributes.get("domain"))
+print("use adaptivity = ",adaptivity_params_ref.attributes.get("use_adaptivity"))
+
+u0_1 = ngs.CoefficientFunction((ngs.x/ngs.sqrt(ngs.x*ngs.x+ngs.y*ngs.y)*ngs.sin(0.5*math.pi*(ngs.x*ngs.x+ngs.y*ngs.y)  )))
+u0_2 = ngs.CoefficientFunction((ngs.y/ngs.sqrt(ngs.x*ngs.x+ngs.y*ngs.y)*ngs.sin(0.5*math.pi*(ngs.x*ngs.x+ngs.y*ngs.y)   )))
+u0_3 = ngs.CoefficientFunction((ngs.cos(0.5*math.pi*(ngs.x*ngs.x+ngs.y*ngs.y) )))
+
+u0 = ngs.CoefficientFunction((u0_1,u0_2,u0_3))
+
+h_mesh = space_init_params_ref.attributes.get("mesh_size")
+print("")
+print("###### tau convergence ######")
+print("mesh size = ", space_init_params_ref.attributes.get("mesh_size"))
+print("Spherical transformation...")
+tau_arr = [1/20,1/40,1/80,1/160,1/320]
+mesh_2D_ref, fes_2D_ref, gfu_2D_ref, gradgfu_2D_ref = SphericalTransform(solver_RSHMHF_ref.gfu[0],solver_RSHMHF_ref.mesh,h_mesh,u0,fes_order)
+print("...done!")
+print("use bdf ", time_init_params_ref.attributes.get("BDF"))
+solver_type = ["PPFEM","TFEM"]
+for s in solver_type:
+    print("")
+    print("\nuse solver ",s)
+    l2errors = []
+    h1errors = []
+    params["time"].attributes["solver_type"] = s
+    for tau in tau_arr:
+        solver_2D = SolvingHMHF2D(u0,params)
+        solver_2D.ChangeTimeStep(tau)
+        print("\ntime step = ",tau)
+        solver_2D.Solve()
+    
+        if s == "TFEM":
+            gradgfu = ngs.CoefficientFunction((ngs.grad(solver_2D.gfu[0].components[0])[0],ngs.grad(solver_2D.gfu[0].components[1])[0],ngs.grad(solver_2D.gfu[0].components[2])[0],ngs.grad(solver_2D.gfu[0].components[0])[1],ngs.grad(solver_2D.gfu[0].components[1])[1],ngs.grad(solver_2D.gfu[0].components[2])[1]), dims=(2,3))
+    
+            l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,ngs.CoefficientFunction((solver_2D.gfu[0].components[0],solver_2D.gfu[0].components[1],solver_2D.gfu[0].components[2])),gradgfu_2D_ref,gradgfu,solver_2D.mesh)
+        else:
+            l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,solver_2D.gfu[0],gradgfu_2D_ref,ngs.grad(solver_2D.gfu[0]),solver_2D.mesh)
+    
+        l2errors.append(l2error)
+        h1errors.append(h1error)
+
+    print("\nSummary")
+    print("time steps = ",tau_arr)
+    eoc_l2 = [ngs.log(l2errors[i+1]/l2errors[i])/ngs.log(0.5) for i in range(len(l2errors)-1)]
+    eoc_h1 = [ngs.log(h1errors[i+1]/h1errors[i])/ngs.log(0.5) for i in range(len(h1errors)-1)]
+    print("L2-errors = ",l2errors)
+    print("L2 eoc = ",eoc_l2)
+    print("")
+    print("H1-errors = ",h1errors)
+    print("H1 eoc = ",eoc_h1)
+
+print("")
+print("###### mesh size convergence ######")
+h_arr = [1/2,1/4,1/8,1/16,1/32,1/64]
+print("time steps = ", time_init_params_ref.attributes.get("tstep"))
+for s in solver_type:
+    print("")
+    print("\nuse solver ",s)
+    params["time"].attributes["solver_type"] = s
+    l2errors = []
+    h1errors = []
+    for h_loop in h_arr:
+    
+        params["space"].attributes["mesh_size"] = h_loop
+        print("\nmesh size = ",h_loop) 
+        print("Spherical transformation...")
+        mesh_2D_ref, fes_2D_ref, gfu_2D_ref, gradgfu_2D_ref = SphericalTransform(solver_RSHMHF_ref.gfu[0],solver_RSHMHF_ref.mesh,h_loop,u0,fes_order)
+        print("...done!")
+        solver_2D = SolvingHMHF2D(u0,params)
+        solver_2D.Solve()
+    
+        if s == "TFEM":
+            gradgfu = ngs.CoefficientFunction((ngs.grad(solver_2D.gfu[0].components[0])[0],ngs.grad(solver_2D.gfu[0].components[1])[0],ngs.grad(solver_2D.gfu[0].components[2])[0],ngs.grad(solver_2D.gfu[0].components[0])[1],ngs.grad(solver_2D.gfu[0].components[1])[1],ngs.grad(solver_2D.gfu[0].components[2])[1]), dims=(2,3))
+    
+            l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,ngs.CoefficientFunction((solver_2D.gfu[0].components[0],solver_2D.gfu[0].components[1],solver_2D.gfu[0].components[2])),gradgfu_2D_ref,gradgfu,solver_2D.mesh)
+        else:
+            l2error,h1error = ComputeH1Error2DForCF(gfu_2D_ref,solver_2D.gfu[0],gradgfu_2D_ref,ngs.grad(solver_2D.gfu[0]),solver_2D.mesh)
+    
+        l2errors.append(l2error)
+        h1errors.append(h1error)
+
+    print("\nSummary")
+    print("mesh sizes = ",h_arr)
+    eoc_l2 = [ngs.log(l2errors[i+1]/l2errors[i])/ngs.log(0.5) for i in range(len(l2errors)-1)]
+    eoc_h1 = [ngs.log(h1errors[i+1]/h1errors[i])/ngs.log(0.5) for i in range(len(h1errors)-1)]
+    print(l2errors)
+    print("L2-errors = ",l2errors)
+    print("L2 eoc = ",eoc_l2)
+    print("")
+    print("H1-errors = ",h1errors)
+    print("H1 eoc = ",eoc_h1)
+
diff --git a/Tests/conv_RSHMHF_FE.py b/Tests/conv_RSHMHF_FE.py
new file mode 100644
index 0000000..fc6fa1e
--- /dev/null
+++ b/Tests/conv_RSHMHF_FE.py
@@ -0,0 +1,158 @@
+## Add project root to sys.path
+import sys
+import os
+sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))
+
+from RSHMHF_FE.solving_RSHMHF import *
+import os
+import numpy as np
+
+## initial condition
+## requirement: u_0(x=0) = 0
+u0 = ngs.CoefficientFunction((math.pi*ngs.x*(1-ngs.x)))
+
+print("############### Compute reference solution ####################")
+
+## input file
+this_dir = os.path.dirname(__file__)
+input_path = os.path.join(this_dir, "input_files/RSHMHF_parameters.txt")
+time_init_params_ref,space_init_params_ref,output_params_ref, adaptivity_params_ref = ReadAndPrintAllParametersFromFile(input_path,0)
+
+number_elements_ref = space_init_params_ref.attributes.get("number_elements")
+use_bdf_ref = time_init_params_ref.attributes.get("BDF")
+dt_ref = time_init_params_ref.attributes.get("tstep")
+fes_order_ref = space_init_params_ref.attributes.get("fes_order")
+
+## ToDo: Put all printing into some outputfile with output options
+print("use bdf ", use_bdf_ref)
+print("time steps = ",dt_ref)
+print("ref number elements = ",number_elements_ref)
+print("fes order = ",fes_order_ref)
+print("use linearization = ", time_init_params_ref.attributes.get("solver_type"))
+print("use adaptivity = ",adaptivity_params_ref.attributes.get("use_adaptivity"))
+a = space_init_params_ref.attributes.get("interval_start")
+b = space_init_params_ref.attributes.get("interval_end")
+mesh_ref = Create1DMesh(number_elements_ref,a,b)
+
+use_dirichlet = space_init_params_ref.attributes.get("Dirichlet")
+if use_dirichlet:
+    ## by default: "left" and "right"
+	dirichlet_str = mesh_ref.GetBoundaries()[0] + "|" + mesh_ref.GetBoundaries()[1]
+	#print(dirichlet_str)
+	V_ref = ngs.H1(mesh_ref,order=fes_order_ref,dirichlet=dirichlet_str)
+	#print(V_ref.FreeDofs())
+else:
+	raise RuntimeError("Only with Dirichlet boundary condition!")
+
+print("ndof = ",V_ref.ndof)
+params_ref = {"time":time_init_params_ref,"output":output_params_ref,"adaptivity":adaptivity_params_ref}
+solver_RSHMHF_ref = SolvingRSHMHF(mesh_ref,V_ref,u0,params_ref)
+solver_RSHMHF_ref.Solve()
+
+print("\n############### Compute solution ####################")
+print("")
+time_init_params,space_init_params,output_params, adaptivity_params = ReadAndPrintAllParametersFromFile(input_path,0)
+
+use_bdf = 1
+time_init_params.attributes["BDF"] = use_bdf
+print("###### Mesh size convergence ######")
+print("use bdf ", use_bdf)
+solver_type = ["convex_lin","max_lin"]
+fes_order = 1
+space_init_params.attributes["fes_order"] = fes_order
+print("\nuse fes order ", fes_order)
+nmb_elements = [2,4,8,16]
+for s in solver_type:
+    l2rerrors = []
+    h1rerrors = []
+    energyerrors = []
+    time_init_params.attributes["solver_type"] = s
+    print("\nuse solver type ",s)
+    for i in nmb_elements:
+    	space_init_params.attributes["number_elements"] = i
+    	use_bdf = time_init_params.attributes.get("BDF")
+    	print("\nmesh size = ",1/i)
+    
+    	mesh = Create1DMesh(i,0,1)
+    	use_dirichlet = space_init_params.attributes.get("Dirichlet")
+    	if use_dirichlet:
+    		dirichlet_str = mesh.GetBoundaries()[0] + "|" + mesh.GetBoundaries()[1]
+    		V = ngs.H1(mesh,order=fes_order,dirichlet=dirichlet_str)
+    		#print(V_ref.FreeDofs())
+    	else:
+    		raise RuntimeError("Only with Dirichlet boundary condition for now!")
+    	
+    	params = {"time":time_init_params,"output":output_params,"adaptivity":adaptivity_params}
+    	solver_RSHMHF = SolvingRSHMHF(mesh,V,u0,params)
+    	solver_RSHMHF.Solve()
+    
+    	l2rerror,h1rerror,energyerror = ComputeH1RAndEnergyError(solver_RSHMHF.gfu[0],solver_RSHMHF_ref.gfu[0],solver_RSHMHF.mesh)
+    	l2rerrors.append(l2rerror)
+    	h1rerrors.append(h1rerror)
+    	energyerrors.append(energyerror)
+    
+    eoc_l2r = [ngs.log(l2rerrors[i+1]/l2rerrors[i])/ngs.log(0.5) for i in range(len(l2rerrors)-1)]
+    eoc_h1r = [ngs.log(h1rerrors[i+1]/h1rerrors[i])/ngs.log(0.5) for i in range(len(h1rerrors)-1)]
+    eoc_energy = [ngs.log(energyerrors[i+1]/energyerrors[i])/ngs.log(0.5) for i in range(len(energyerrors)-1)]
+    print("\nSummary")
+    print("mesh sizes = ",[1/2,1/4,1/8,1/16])
+    print("L2-errors = ",l2rerrors)
+    print("L2 eoc = ",eoc_l2r)
+    print("")
+    print("H1-errors = ",h1rerrors)
+    print("H1 eoc",eoc_h1r)
+    print("")
+    print("Energy-errors = ",energyerrors)
+    print("Energy eoc = ",eoc_energy)
+
+print("\n###### Time step convergence ######")
+print("")
+time_init_params.attributes["BDF"] = 1
+space_init_params.attributes["number_elements"] =number_elements_ref 
+dt_arr = [1.25e-2,6.25e-3,3.125e-3,3.125e-3/2,3.125e-3/4]
+
+space_init_params.attributes["fes_order"] = fes_order
+print("use fes order = ",fes_order)
+print("use bdf ", use_bdf)
+print("use mesh size = ",1/number_elements_ref)
+mesh = Create1DMesh(number_elements_ref,0,1)
+use_dirichlet = space_init_params.attributes.get("Dirichlet")
+if use_dirichlet:
+    dirichlet_str = mesh.GetBoundaries()[0] + "|" + mesh.GetBoundaries()[1]
+    V = ngs.H1(mesh,order=fes_order,dirichlet=dirichlet_str)
+	#print(V_ref.FreeDofs())
+else:
+    raise RuntimeError("Only with Dirichlet boundary condition for now!")
+for s in solver_type:
+    l2rerrors = []
+    h1rerrors = []
+    energyerrors = []
+    time_init_params.attributes["solver_type"] = s
+    print("\nuse solver type ",s)
+    for dt in dt_arr:
+        time_init_params.attributes["tstep"] = dt
+        time_step = time_init_params.attributes.get("tstep")
+        print("\ntime steps = ",time_step)
+    	
+        params = {"time":time_init_params,"output":output_params,"adaptivity":adaptivity_params}
+        solver_RSHMHF = SolvingRSHMHF(mesh,V,u0,params)
+        solver_RSHMHF.Solve()
+    
+        l2rerror,h1rerror,energyerror = ComputeH1RAndEnergyError(solver_RSHMHF.gfu[0],solver_RSHMHF_ref.gfu[0],solver_RSHMHF.mesh)
+        l2rerrors.append(l2rerror)
+        h1rerrors.append(h1rerror)
+        energyerrors.append(energyerror)
+    
+    eoc_l2r = [ngs.log(l2rerrors[i+1]/l2rerrors[i])/ngs.log(0.5) for i in range(len(l2rerrors)-1)]
+    eoc_h1r = [ngs.log(h1rerrors[i+1]/h1rerrors[i])/ngs.log(0.5) for i in range(len(h1rerrors)-1)]
+    eoc_energy = [ngs.log(energyerrors[i+1]/energyerrors[i])/ngs.log(0.5) for i in range(len(energyerrors)-1)]
+    print("\nSummary")
+    print("time steps = ",dt_arr)
+    print("L2-errors = ",l2rerrors)
+    print("L2 eoc = ",eoc_l2r)
+    print("")
+    print("H1-errors = ",h1rerrors)
+    print("H1 eoc = ",eoc_h1r)
+    print("")
+    print("Energy-errors = ",energyerrors)
+    print("Energy eoc = ",eoc_energy)
diff --git a/Tests/input_files/HMHF_2D_PPFEM_parameters.txt b/Tests/input_files/HMHF_2D_PPFEM_parameters.txt
new file mode 100644
index 0000000..b398c53
--- /dev/null
+++ b/Tests/input_files/HMHF_2D_PPFEM_parameters.txt
@@ -0,0 +1,28 @@
+# for now only unit disk allowed
+domain unit_disk
+
+# number of elements in interval I
+mesh_size 0.00390625
+ 
+# time parameters
+tstart 0.0
+tstep 1e-6
+tend 0.1
+
+# time stepping scheme
+BDF 1
+
+# 1 = true, 0 = false
+# if 0, then zero Neumann used
+Dirichlet 1
+
+Dirichlet_bnd_name unit_circle
+
+fes_order 1
+
+# debug mode: print to terminal
+debug 1
+
+# PPFEM, TFEM or CPFEM
+# CPFEM only convergence in mesh size feasible
+solver_type PPFEM
diff --git a/Tests/input_files/RSHMHF_parameters.txt b/Tests/input_files/RSHMHF_parameters.txt
new file mode 100644
index 0000000..d89f236
--- /dev/null
+++ b/Tests/input_files/RSHMHF_parameters.txt
@@ -0,0 +1,28 @@
+# define the intervall in 1D: I := [interval_start,interval_end]
+interval_start 0
+interval_end 1
+
+# number of elements in interval I, h = 1/number_elements
+number_elements 16384
+#number_elements 10
+
+# different linearizations for non-linearity in RSHMHF 
+solver_type convex_lin
+
+# time parameters
+tstart 0.0
+tstep 1e-6
+#tstep 1e-2
+tend 0.1
+
+# time stepping scheme
+BDF 2
+
+fes_order 2
+
+# 1 = true, 0 = false
+# if 0, then zero Neumann used
+Dirichlet 1
+
+# debug mode: print to terminal
+debug 1
-- 
GitLab