--- a/src/SbpOperators/boundary_conditions/boundary_condition.jl	Sun Jun 09 17:04:01 2024 -0700
+++ b/src/SbpOperators/boundary_conditions/boundary_condition.jl	Sun Jun 09 17:51:57 2024 -0700
@@ -1,7 +1,7 @@
 """
     BoundaryCondition{BID}
 
-A type for implementing boundary_data needed in order to impose a boundary condition.
+A type for implementing data needed in order to impose a boundary condition.
 Subtypes refer to perticular types of boundary conditions, e.g. Neumann conditions.
 """
 abstract type BoundaryCondition{BID} end
@@ -30,6 +30,12 @@
     return eval_on(boundary_grid(grid, boundary(bc)), boundary_data(bc))
 end
 
+"""
+    DirichletCondition{DT,BID}
+
+A Dirichlet condition with `data::DT` on the boundary
+specified by the boundary identifier `BID`.
+"""
 struct DirichletCondition{DT,BID} <: BoundaryCondition{BID}
     data::DT
     function DirichletCondition(data, id)
@@ -38,6 +44,12 @@
 end
 boundary_data(bc::DirichletCondition) = bc.data
 
+"""
+    DirichletCondition{DT,BID}
+
+A Neumann condition with `data::DT` on the boundary
+specified by the boundary identifier `BID`.
+"""
 struct NeumannCondition{DT,BID} <: BoundaryCondition{BID}
     data::DT
     function NeumannCondition(data, id)

--- a/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 09 17:04:01 2024 -0700
+++ b/src/SbpOperators/volumeops/laplace/laplace.jl	Sun Jun 09 17:51:57 2024 -0700
@@ -53,14 +53,13 @@
 end
 laplace(g::EquidistantGrid, stencil_set) = second_derivative(g, stencil_set)
 
-
 """
-sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; tuning)
+    sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
 
 The operators required to construct the SAT for imposing a Dirichlet condition.
-`tuning` specifies the strength of the penalty. See
+`H_tuning` and `R_tuning` are used to specify the strength of the penalty.
 
-See also: [`sat`,`DirichletCondition`, `positivity_decomposition`](@ref).
+See also: [`sat`](@ref),[`DirichletCondition`](@ref),[`positivity_decomposition`](@ref).
 """
 function sat_tensors(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning = 1., R_tuning = 1.)
     id = boundary(bc)
@@ -75,12 +74,11 @@
 end
 
 """
-sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
+    sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
 
 The operators required to construct the SAT for imposing a Neumann condition
 
-
-See also: [`sat`,`NeumannCondition`](@ref).
+See also: [`sat`](@ref),[`NeumannCondition`](@ref).
 """
 function sat_tensors(Δ::Laplace, g::Grid, bc::NeumannCondition)
     id = boundary(bc)
@@ -94,13 +92,22 @@
     return penalty_tensor, d
 end
 
+"""
+    positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
 
+Constructs the scalar `B` such that `d' - 1/2*B*e'` is symmetric positive definite with respect to
+the boundary quadrature. Here `d` is the normal derivative and `e` is the boundary restriction operator.
+`B` can then be used to form a symmetric and energy stable penalty for a Dirichlet condition.
+The parameters `H_tuning` and `R_tuning` are used to specify the strength of the penalty and
+must be greater than 1. For details we refer to https://doi.org/10.1016/j.jcp.2020.109294
+"""
 function positivity_decomposition(Δ::Laplace, g::Grid, bc::DirichletCondition; H_tuning, R_tuning)
+    @assert(H_tuning ≥ 1.)
+    @assert(R_tuning ≥ 1.)
     Nτ_H, τ_R = positivity_limits(Δ,g,bc)
     return H_tuning*Nτ_H + R_tuning*τ_R
 end
 
-
 # TODO: We should consider implementing a proper BoundaryIdentifier for EquidistantGrid and then
 # change bc::BoundaryCondition to id::BoundaryIdentifier
 function positivity_limits(Δ::Laplace, g::EquidistantGrid, bc::DirichletCondition)