Repository: magpylib/magpylib (language: python)

Issue:
The Triangle class
Hi everyone, but specifically @Alexboiboi @fslanovc  and @peter-leitner 

After a lot of thinking and discussion with you I propose the following interface layout for the Facet code: We introduce one core function `magnet_facet_field` and three classes for the object oriented interface: `magnet.Tetrahedron` , `magnet.FacetBody`, and `misc.Facet`

The Facet class
-------------------
Parameters: `magnetization` (the usual) and `vertices` shape (3,3) or (n,3,3)

This class contains only the bare triangular facet code and is for users who like to customize, e.g. work only with single faces, or non-closed meshes because they know that parallel faces to the magnetization do not contribute.

Users have to think by themselves about inside/outside, and about the face orientations.

Facet objects are not magnets (not even if all faces are properly arranged, because there is no inside-outside logic), so the `Facet` goes to the `misc` sub-package.

This already works in the [facet_MO](https://github.com/magpylib/magpylib/tree/facet_MO) branch 

The Tetrahedron class
-------------------------
Parameters: `magnetization` (the usual) and `vertices` shape (4,3)

This class is introduced because the Tetrahedron is a very basic shape which can be used to construct a lot of complex forms.

The Tetrahedron is very user-friendly: The faces are automatically arranged in the computation and inside-outside is computed. It behaves like a magnet, so it goes to the `magnet` sub-package.

This already works in the [facet_MO](https://github.com/magpylib/magpylib/tree/facet_MO) branch 



Discussion / hints:

Fix patch:
diff --git a/.gitignore b/.gitignore
index 810fdda59..329990374 100644
--- a/.gitignore
+++ b/.gitignore
@@ -37,6 +37,8 @@ __temp.py
 __temp2.py
 __temp3.py
 __temp4.py
+__temp5.py
+__temp6.py
 
 # PyInstaller
 #  Usually these files are written by a python script from a template
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2a7ff5652..1cbab8e24 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,6 +2,19 @@ All notable changes to magpylib are documented here.
 
 
 # Releases
+
+## [4.2.0] - XXXXXXXXXXX
+- (Re)introducing the powerful `misc.Triangle` class that can be used to compute magnetic fields of arbitrarily shaped bodies by approximating their surface with triangular facets. ([#568](https://github.com/magpylib/magpylib/issues/568))
+- Introducing the `magnet.Tetrahedron` class as a derivate of the Triangle class. ([#289](https://github.com/magpylib/magpylib/issues/289))
+- Added code of conduct attempting to align with NumFocus standards ([#558](https://github.com/magpylib/magpylib/issues/558))
+- Improved Loop field computation in terms of performance and numerical stability ([#374](https://github.com/magpylib/magpylib/issues/374))
+- Documentation changes:
+    - Correct conda install command
+    - Integration of Triangle and Tetrahedron
+    - Changed example gallery substructure
+    - Rewritten some passages
+- Fixed some bugs, minor performance increase, internal refactoring
+
 ## [4.1.2] - 2023-01-15
 - Fix wrong magnetization arrow direction for some edge cases ([#570](https://github.com/magpylib/magpylib/discussions/570), [#571](https://github.com/magpylib/magpylib/issues/571), [#572](https://github.com/magpylib/magpylib/pull/572))
 - Fix cryptic `getB`/`getH` error message ([#562](https://github.com/magpylib/magpylib/issues/562), [#563](https://github.com/magpylib/magpylib/pull/563))
diff --git a/CODE_OF_CONDUCT.md b/CODE_OF_CONDUCT.md
new file mode 100644
index 000000000..e0c7ab667
--- /dev/null
+++ b/CODE_OF_CONDUCT.md
@@ -0,0 +1,73 @@
+# Magpylib Code of Conduct
+
+## Our Pledge
+
+We as members, contributors, and leaders pledge to make participation in our community a harassment-free experience for everyone, regardless of age, body size, visible or invisible disability, ethnicity, sex characteristics, gender identity and expression, level of experience, education, socio-economic status, nationality, personal appearance, race, caste, color, religion, or sexual identity and orientation.
+
+We pledge to act and interact in ways that contribute to an open, welcoming, diverse, inclusive, and healthy community.
+
+## Our Standards
+
+Examples of behavior that contributes to a positive environment for our community include:
+
+* Demonstrating empathy and kindness toward other people
+* Being respectful of differing opinions, viewpoints, and experiences
+* Giving and gracefully accepting constructive feedback
+* Accepting responsibility and apologizing to those affected by our mistakes, and learning from the experience
+* Focusing on what is best not just for us as individuals, but for the overall community
+
+Examples of unacceptable behavior include:
+
+* The use of sexualized language or imagery, and sexual attention or advances of any kind
+* Trolling, insulting or derogatory comments, and personal or political attacks
+* Public or private harassment
+* Publishing others' private information, such as a physical or email address, without their explicit permission
+* Other conduct which could reasonably be considered inappropriate in a professional setting
+
+## Enforcement Responsibilities
+
+Community leaders are responsible for clarifying and enforcing our standards of acceptable behavior and will take appropriate and fair corrective action in response to any behavior that they deem inappropriate, threatening, offensive, or harmful.
+
+Community leaders have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct, and will communicate reasons for moderation decisions when appropriate.
+
+## Scope
+
+This Code of Conduct applies within all community spaces, and also applies when an individual is officially representing the community in public spaces. Examples of representing our community include using an official e-mail address, posting via an official social media account, or acting as an appointed representative at an online or offline event.
+
+## Enforcement
+
+Instances of abusive, harassing, or otherwise unacceptable behavior may be reported to the community leaders responsible for enforcement at [magpylib@gmail.com](mailto:magpylib@gmail.com). All complaints will be reviewed and investigated promptly and fairly.
+
+All community leaders are obligated to respect the privacy and security of the reporter of any incident.
+
+## Enforcement Guidelines
+
+Community leaders will follow these Community Impact Guidelines in determining the consequences for any action they deem in violation of this Code of Conduct:
+
+### 1. Correction
+
+**Community Impact**: Use of inappropriate language or other behavior deemed unprofessional or unwelcome in the community.
+
+**Consequence**: A private, written warning from community leaders, providing clarity around the nature of the violation and an explanation of why the behavior was inappropriate. A public apology may be requested.
+
+### 2. Warning
+
+**Community Impact**: A violation through a single incident or series of actions.
+
+**Consequence**: A warning with consequences for continued behavior. No interaction with the people involved, including unsolicited interaction with those enforcing the Code of Conduct, for a specified period of time. This includes avoiding interactions in community spaces as well as external channels like social media. Violating these terms may lead to a temporary or permanent ban.
+
+### 3. Temporary Ban
+
+**Community Impact**: A serious violation of community standards, including sustained inappropriate behavior.
+
+**Consequence**: A temporary ban from any sort of interaction or public communication with the community for a specified period of time. No public or private interaction with the people involved, including unsolicited interaction with those enforcing the Code of Conduct, is allowed during this period. Violating these terms may lead to a permanent ban.
+
+### 4. Permanent Ban
+
+**Community Impact**: Demonstrating a pattern of violation of community standards, including sustained inappropriate behavior, harassment of an individual, or aggression toward or disparagement of classes of individuals.
+
+**Consequence**: A permanent ban from any sort of public interaction within the community.
+
+## Attribution
+
+This Code of Conduct is adapted from the [Contributor Covenant](https://www.contributor-covenant.org/).
\ No newline at end of file
diff --git a/docs/_pages/page_01_introduction.md b/docs/_pages/page_01_introduction.md
index 568deff39..ce01f03ea 100644
--- a/docs/_pages/page_01_introduction.md
+++ b/docs/_pages/page_01_introduction.md
@@ -15,7 +15,7 @@ kernelspec:
 
 # Introduction
 
-This section gives an introduction to the Magpylib API. More detailed information and practical examples how to use Magpylib can be found in the example galleries. The package documentation is found in the {ref}`modindex`.
+This section gives an introduction to the Magpylib API. More detailed information and practical examples how to use Magpylib can be found in the example galleries. The package docstrings are found in the {ref}`modindex`.
 
 ## Contents
 
@@ -33,10 +33,12 @@ This section gives an introduction to the Magpylib API. More detailed informatio
 
 ## The idea behind Magpylib
 
-Magpylib provides fast and accurate magnetostatic field computation based on **analytical solutions** to permanent magnet and current problems. The field computation is coupled to a **position and orientation interface** that makes it easy to work with relative object positioning. The idea behind the main object oriented interface is:
+Magpylib provides fast and accurate magnetostatic field computation based on **analytical solutions** to permanent magnet and current problems. The fast field computation is coupled to a **position and orientation API** that makes it easy to work with relative object positioning.
+
+The idea behind the main object oriented interface is:
 
 1. Sensors, magnets, currents, etc. are created as Python objects with defined position and orientation in a global coordinate system.
-2. After initialization, the Magpylib objects can easily be manipulated, grouped, moved around and displayed graphically using Matplotlib or Plotly 3D plotting.
+2. After initialization, the Magpylib objects can easily be manipulated, grouped, moved around and displayed graphically using 3D graphic backends (Matplotlib, Plotly, Pyvista).
 3. When all objects are set up, the magnetic field generated by the source objects is computed at the observer objects.
 
 The following example code outlines this functionality:
@@ -75,7 +77,7 @@ For users who would like to avoid the object oriented interface, the field imple
 
 The analytical solutions are exact when there is no material response and natural boundary conditions can be assumed. In general, Magpylib is at its best when dealing with air-coils (no eddy currents) and high grade permanent magnets (Ferrite, NdFeB, SmCo or similar materials).
 
-When **magnet** permeabilities are below $\mu_r < 1.1$ the error typically undercuts 1-5 % (long magnet shapes are better, large distance from magnet is better). Demagnetization factors are not automatically included at this point. The line **current** solutions give the exact same field as outside of a wire that carries a homogenous current. For more details check out the {ref}`physComp` section.
+When **magnet** permeabilities are below $\mu_r < 1.1$ the error typically undercuts 1-5 % (long magnet shapes are better, large distance from magnet is better). Demagnetization factors are not automatically included at this point. They typically reduce the error by an order of magnitude. The line **current** solutions give the exact same field as outside of a wire that carries a homogenous current. For more details check out the {ref}`physComp` section.
 
 Magpylib only provides solutions for simple geometric forms (cuboids, cylinders, lines, ...). How **complex shapes** can be constructed from these simple base shapes is described in {ref}`examples-complex-forms`.
 
@@ -89,53 +91,61 @@ The most convenient way of working with Magpylib is through the **object oriente
 
 All magnet objects have the `magnetization` attribute which must be of the format $(m_x, m_y, m_z)$ and denotes the homogeneous magnetization vector in the local object coordinates in units of \[mT\]. It is often referred to as the remanence ($B_r=\mu_0 M$) in material data sheets. All magnets can be used as Magpylib `sources` input.
 
-- **`Cuboid`**`(magnetization, dimension, position, orientation)` represents a magnet with cuboid shape. `dimension` has the format $(a,b,c)$ and denotes the sides of the cuboid in units of \[mm\]. By default the center of the cuboid lies in the origin of the global coordinates, and the sides are parallel to the coordinate axes.
+- **`Cuboid`**`(magnetization, dimension, position, orientation, style)` represents a magnet with cuboid shape. `dimension` has the format $(a,b,c)$ and denotes the sides of the cuboid in units of \[mm\]. By default the center of the cuboid lies in the origin of the global coordinates, and the sides are parallel to the coordinate axes.
+
+- **`Cylinder`**`(magnetization, dimension, position, orientation, style)` represents a magnet with cylindrical shape. `dimension` has the format $(d,h)$ and denotes diameter and height of the cylinder in units of \[mm\]. By default the center of the cylinder lies in the origin of the global coordinates, and the cylinder axis coincides with the z-axis.
 
-- **`Cylinder`**`(magnetization, dimension, position, orientation)` represents a magnet with cylindrical shape. `dimension` has the format $(d,h)$ and denotes diameter and height of the cylinder in units of \[mm\]. By default the center of the cylinder lies in the origin of the global coordinates, and the cylinder axis coincides with the z-axis.
+- **`CylinderSegment`**`(magnetization, dimension, position, orientation, style)` represents a magnet with the shape of a cylindrical ring section. `dimension` has the format $(r_1,r_2,h,\varphi_1,\varphi_2)$ and denotes inner radius, outer radius and height in units of \[mm\] and the two section angles $\varphi_1<\varphi_2$ in \[deg\]. By default the center of the full cylinder lies in the origin of the global coordinates, and the cylinder axis coincides with the z-axis.
 
-- **`CylinderSegment`**`(magnetization, dimension, position, orientation)` represents a magnet with the shape of a cylindrical ring section. `dimension` has the format $(r_1,r_2,h,\varphi_1,\varphi_2)$ and denotes inner radius, outer radius and height in units of \[mm\] and the two section angles $\varphi_1<\varphi_2$ in \[deg\]. By default the center of the full cylinder lies in the origin of the global coordinates, and the cylinder axis coincides with the z-axis.
+- **`Sphere`**`(magnetization, diameter, position, orientation, style)` represents a magnet of spherical shape. `diameter` is the sphere diameter $d$ in units of \[mm\]. By default the center of the sphere lies in the origin of the global coordinates.
 
-- **`Sphere`**`(magnetization, diameter, position, orientation)` represents a magnet of spherical shape. `diameter` is the sphere diameter $d$ in units of \[mm\]. By default the center of the sphere lies in the origin of the global coordinates.
+- **`Tetrahedron`**`(magnetization, vertices, position, orientation, style)` represents a magnet of tetrahedral shape. `vertices` corresponds to the four corner points in units of \[mm\]. By default the vertex positions coincide in the local object coordinates and the global coordinates.
 
 **Currents**
 
 All current objects have the `current` attribute which must be a scalar $i_0$ and denotes the electrical current in units of \[A\]. All currents can be used as Magpylib `sources` input.
 
-- **`Loop`**`(current, diameter, position, orientation)` represents a circular current loop where `diameter` is the loop diameter $d$ in units of \[mm\]. By default the loop lies in the xy-plane with it's center in the origin of the global coordinates.
+- **`Loop`**`(current, diameter, position, orientation, style)` represents a circular current loop where `diameter` is the loop diameter $d$ in units of \[mm\]. By default the loop lies in the xy-plane with it's center in the origin of the global coordinates.
 
-- **`Line`**`(current, vertices, position, orientation)` represents electrical current segments that flow in a straight line from vertex to vertex. By default the verticx positions coincide in the local object coordinates and the global coordinates.
+- **`Line`**`(current, vertices, position, orientation, style)` represents electrical current segments that flow in a straight line from vertex to vertex. By default the vertex positions coincide in the local object coordinates and the global coordinates.
 
 **Other**
 
-- **`Dipole`**`(moment, position, orientation)` represents a magnetic dipole moment with moment $(m_x,m_y,m_z)$ given in \[mT mm³]. For homogeneous magnets the relation moment=magnetization$\times$volume holds. Can be used as Magpylib `sources` input.
+- **`Dipole`**`(moment, position, orientation, style)` represents a magnetic dipole moment with moment $(m_x,m_y,m_z)$ given in \[mT mm³]. For homogeneous magnets the relation moment=magnetization$\times$volume holds. Can be used as Magpylib `sources` input.
+
+- **`Triangle`**`(magnetization, vertices, position, orientation, style)` represents a triangular surface with a homogeneous charge given by the projection of the `magnetization` vector onto the surface normal. `vertices` is a set of the three corners given in \[mm³\]. When multiple Triangles form a closed surface, on the outside their total magnetic field correponds to the one of a homogeneously charged magnet.
 
-- **`CustomSource`**`(field_func, position, orientation)` is used to create user defined custom sources with their own field functions. Can be used as Magpylib `sources` input.
+- **`CustomSource`**`(field_func, position, orientation, style)` is used to create user defined custom sources with their own field functions. Can be used as Magpylib `sources` input.
 
-- **`Sensor`**`(position, pixel, orientation)` represents a 3D magnetic field sensor. The field is evaluated at the given pixel positions. By default (`pixel=(0,0,0)`) the pixel position coincide in the local object coordinates and the global coordinates. Can be used as Magpylib `observers` input.
+- **`Sensor`**`(position, pixel, orientation, style)` represents a 3D magnetic field sensor. The field is evaluated at the given pixel positions. By default (`pixel=(0,0,0)`) the pixel position coincide in the local object coordinates and the global coordinates. Can be used as Magpylib `observers` input.
 
-- **`Collection`**`(*children, position, orientation)` is a group of source and sensor objects (children) that is used for common manipulation. Depending on the children, a collection can be used as Magpylib `sources` and/or `observers` input.
+- **`Collection`**`(*children, position, orientation, style)` is a group of source and sensor objects (children) that is used for common manipulation. Depending on the children, a collection can be used as Magpylib `sources` and/or `observers` input.
 
 ```{code-cell} ipython3
 import magpylib as magpy
 
+objects = []
+
 # magnets
-magnet1 = magpy.magnet.Cuboid()
-magnet2 = magpy.magnet.Cylinder()
-magnet3 = magpy.magnet.CylinderSegment()
-magnet4 = magpy.magnet.Sphere()
+objects.append(magpy.magnet.Cuboid())
+objects.append(magpy.magnet.Cylinder())
+objects.append(magpy.magnet.CylinderSegment())
+objects.append(magpy.magnet.Sphere())
+objects.append(magpy.magnet.Tetrahedron())
 
 # currents
-current1 = magpy.current.Loop()
-current2 = magpy.current.Line()
+objects.append(magpy.current.Loop())
+objects.append(magpy.current.Line())
 
 # other
-dipole = magpy.misc.Dipole()
-custom = magpy.misc.CustomSource()
-sensor = magpy.Sensor()
-coll   = magpy.Collection()
+objects.append(magpy.misc.Dipole())
+objects.append(magpy.misc.Triangle())
+objects.append(magpy.misc.CustomSource())
+objects.append(magpy.Sensor())
+objects.append(magpy.Collection())
 
 # print object representation
-for obj in [magnet1, magnet2, magnet3, magnet4, current1, current2, dipole, custom, sensor, coll]:
+for obj in objects:
     print(obj)
 ```
 
@@ -143,24 +153,24 @@ for obj in [magnet1, magnet2, magnet3, magnet4, current1, current2, dipole, cust
 
 ## Position and orientation
 
-All Magpylib objects have the `position` and `orientation` attributes that refer to position and orientation in the global Cartesian coordinate system. The `position` attribute is a numpy ndarray, shape (3,) or (m,3) and denotes the coordinates $(x,y,z)$ in units of [mm]. By default every object is created at `position=(0,0,0)`. The `orientation` attribute is a scipy [Rotation object](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.html) and denotes the object rotation relative to its initial state. By default the orientation of a Magpylib object is the unit rotation, `orientation=None`.
+All Magpylib objects have the `position` and `orientation` attributes that refer to position and orientation in the global Cartesian coordinate system. The `position` attribute is a numpy ndarray, shape (3,) or (m,3) and corresponds to the coordinates $(x,y,z)$ in units of [mm]. By default every object is created in the origin of the global coordinate system, at `position=(0,0,0)`. The `orientation` attribute is a scipy [Rotation object](https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.Rotation.html) and corresponds to the object orientation relative to its initial state. By default objects are created with initial orientation, `orientation=None`, which is the unit rotation. The initial orientation of every object, e.g. current loop lies in the x-y plane, is described in the respective class docstrings.
 
 ```python
 import magpylib as magpy
 
+# init object with default values
 sensor = magpy.Sensor()
-
 print(sensor.position)                                     # out: [0. 0. 0.]
 print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [0. 0. 0.]
 ```
 
-Set absolute object position/orientation at initialization or through the properties directly. Add relative position/orientation with the `move` and `rotate` methods,
+Set absolute object position and orientation attributes at object initialization or directly through the properties.
 
 ```python
 import magpylib as magpy
 from scipy.spatial.transform import Rotation as R
 
-# set at initialization
+# set attributes at initialization
 sensor = magpy.Sensor(position=(1,1,1))
 print(sensor.position)                                     # out: [1. 1. 1.]
 print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [0. 0. 0.]
@@ -169,12 +179,30 @@ print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [0. 0. 0.]
 sensor.orientation = R.from_rotvec((0,0,45), degrees=True)
 print(sensor.position)                                     # out: [1. 1. 1.]
 print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [ 0.  0. 45.]
+```
+
+The `move` and `rotate` methods are powerful tools to change the relative position and orientation of an existing object. **Move** the object by a `displacement` vector. **Rotate** an object by specifying the angle of rotation `angle` (scalar), an axis of rotation `axis` (vector or `'x'`, `'y'`, `'z'`) and an anchor point `anchor` (vector) through which the rotation axis passes through. By default `anchor=self.position`, meaning that the object rotates about itself.
+
+```python
+import magpylib as magpy
+
+# init object with default values
+sensor = magpy.Sensor()
+
+# move
+sensor.move(displacement=(1,1,3))
+print(sensor.position)                                     # out: [1. 1. 3.]
+print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [ 0.  0.  0.]
+
+# rotate about self
+sensor.rotate_from_angax(angle=45, axis='z')
+print(sensor.position)                                     # out: [1. 1. 3.]
+print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [ 0.  0. 45.]
 
-# move and rotate
-sensor.move((1,2,3))
-sensor.rotate_from_angax(45, 'z')
-print(sensor.position)                                     # out: [2. 3. 4.]
-print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [ 0.  0. 90.]
+# rotate with anchor
+sensor.rotate_from_angax(angle=90, axis='z', anchor=(0,0,0))
+print(sensor.position)                                     # out: [-1. 1. 3.]
+print(sensor.orientation.as_euler('xyz', degrees=True))    # out: [ 0.  0. 135.]
 ```
 
 The attributes `position` and `orientation` can be either of **"scalar"** nature, i.e. a single position or a single rotation like in the examples above, or **"vectors"** when they are arrays of such scalars. The two attributes together define an object **"path"**. Paths should always be used when modelling object motion as the magnetic field is computed on the whole path with increased performance.
@@ -219,9 +247,9 @@ The following example shows the graphical representation of various Magpylib obj
 ```{code-cell} ipython3
 import numpy as np
 import magpylib as magpy
-from magpylib.magnet import Cuboid, Cylinder, CylinderSegment, Sphere
+from magpylib.magnet import Cuboid, Cylinder, CylinderSegment, Sphere, Tetrahedron
 from magpylib.current import Loop, Line
-from magpylib.misc import Dipole
+from magpylib.misc import Dipole, Triangle
 
 objects = [
     Cuboid(
@@ -244,6 +272,11 @@ objects = [
         diameter=1,
         position=(-1,0,0),
     ),
+    Tetrahedron(
+        magnetization=(0,0,100),
+        vertices=((-1,0,0), (1,0,0), (0,-1,0), (0,0,1)),
+        position=(-1,0,3)
+    ),
     Loop(
         current=1,
         diameter=1,
@@ -258,6 +291,11 @@ objects = [
         moment=(0,0,100),
         position=(5,0,0),
     ),
+    Triangle(
+        magnetization=(0,0,100),
+        vertices=((-1,0,0), (1,0,0), (0,1,0)),
+        position=(5,0,3),
+    ),
     magpy.Sensor(
         pixel=[(0,0,z) for z in (-.5,0,.5)],
         position=(7,0,0),
@@ -568,7 +606,7 @@ magpy.show(*helmholtz, backend='plotly', animation=4, style_path_show=False)
 
 Notice, that collections have their own `style` attributes, their paths are displayed in `show`, and all children are automatically assigned their parent color.
 
-For magnetic field computation a collection with source children behaves like a single source object, and a collection with sensor children behaves like a flat list of it's sensors when provided as `sources` and `observers` input respectively. This is demonstrated in the following continuation of the previous helmholtz example:
+For magnetic field computation a collection with source children behaves like a single source object, and a collection with sensor children behaves like a flat list of it's sensors when provided as `sources` and `observers` input respectively. This is demonstrated in the following continuation of the previous Helmholtz example:
 
 ```{code-cell} ipython3
 import matplotlib.pyplot as plt
@@ -616,7 +654,6 @@ egg = magpy.misc.CustomSource(
 
 # add a custom 3D model
 trace = magpy.graphics.model3d.make_Ellipsoid(
-    backend='plotly',
     dimension=(1,1,1.4),
 )
 egg.style.model3d.add_trace(trace)
diff --git a/docs/_pages/page_03_installation.md b/docs/_pages/page_03_installation.md
index ca60179bc..9c315615a 100644
--- a/docs/_pages/page_03_installation.md
+++ b/docs/_pages/page_03_installation.md
@@ -11,7 +11,7 @@ pip install magpylib
 or with [conda](https://docs.conda.io/en/latest/).
 
 ```console
-conda install magpylib
+conda install -c conda-forge magpylib
 ```
 
 Magpylib is hosted for **download** at the following locations:
diff --git a/docs/_pages/page_07_physics_computation.md b/docs/_pages/page_07_physics_computation.md
index c222741cd..de41b12b1 100644
--- a/docs/_pages/page_07_physics_computation.md
+++ b/docs/_pages/page_07_physics_computation.md
@@ -6,15 +6,15 @@
 
 **Permanent Magnets**
 
-Magnetic field computations in Magpylib are based on known analytical solutions (formulas) to permanent magnet and current problems. For Magpylib we have used the following references:
+Magnetic field computations in Magpylib are based on known analytical solutions (closed form expressions) to permanent magnet and current problems. The Magpylib implementations are based on the following literature references:
 
 - Field of cuboid magnets: \[Yang1999, Engel-Herbert2005, Camacho2013, Cichon2019\]
 - Field of cylindrical magnets: \[Furlani1994, Derby2009, Caciagli2018, Slanovc2021\]
-- Field of facet bodies: \[Janssen2009, Rubeck2013\]
-- Field of circular line current: \[Smythe1950, Simpson2001, Ortner2017\]
+- Field of triangular surfaces: \[Guptasarma1999, Janssen2009, Rubeck2013\]
+- Field of the current loop: \[Ortner2022\]
 - all others derived by hand
 
-A short reflection on how these formulas can be achieved: In magnetostatics (no currents) the magnetic field becomes conservative (Maxwell: $\nabla \times {\bf H} = 0$) and can thus be expressed through the magnetic scalar potential $\Phi_m$:
+A short reflection on how these formulas can be achieved: In magnetostatics the magnetic field becomes conservative (Maxwell: $\nabla \times {\bf H} = 0$) and can thus be expressed through the magnetic scalar potential $\Phi_m$:
 
 $$
 {\bf H} = -\nabla\cdot\Phi_m
@@ -26,7 +26,7 @@ $$
 \Phi_m({\bf r}) = \frac{1}{4\pi}\int_{V'}\frac{\nabla'\cdot {\bf M}({\bf r}')}{|{\bf r}-{\bf r}'|}dV'+\frac{1}{4\pi}\oint_{S'}\frac{{\bf n}'\cdot {\bf M}({\bf r}')}{|{\bf r}-{\bf r}'|}dS'
 $$
 
-where ${\bf r}$ denotes the position, $V$ is the magnetized volume with surface $S$ and normal vector ${\bf n}$ onto the surface. This solution is derived in detail e.g. in \[1999Jackson\].
+where ${\bf r}$ denotes the position, $V$ is the magnetized volume with surface $S$ and normal vector ${\bf n}$ onto the surface. This solution is derived in detail e.g. in \[Jackson1999\].
 
 **Currents**
 
@@ -42,7 +42,7 @@ In some special cases (simple shapes, homogeneous magnetizations and current dis
 
 **Line currents:**
 
-The magnetic field of a wire carrying a homogeneous current density is similar (ON THE OUTSIDE ONLY) to the field of a line current in the center of the wire, which carries the total current of the wire. Current distributions become inhomogeneous at bends of the wire or when eddy currents (finite frequencies) are involved.
+The magnetic field of a wire carrying a homogeneous current density is (on the outside) similar to the field of a line current positioned in the center of the wire, which carries the same total current. Current distributions become inhomogeneous at bends of the wire or when eddy currents (finite frequencies) are involved.
 
 **Magnets and Demagnetization**
 
@@ -56,7 +56,7 @@ Demagnetization factors can be used to compensate a large part of the demagnetiz
 
 **Soft-Magnetic Materials**
 
-Soft-magnetic materials like iron or steel with large permeabilities $\mu_r \sim 1000$ and low remanence fields are dominated by the material response. There are no useful analytical solutions. However, recent developments showed that the Magnetostatic Method of Moments can be a powerful tool in combination with Magpylib to compute such a material response. An integration into Magpylib is planned for the future.
+Soft-magnetic materials like iron or steel with large permeabilities $\mu_r \sim 1000$ and low remanence fields are dominated by the material response. It is not possible to describe such bodies with analytical solutions. However, recent developments showed that the Magnetostatic Method of Moments can be a powerful tool in combination with Magpylib to compute such a material response. An integration into Magpylib is planned for the future.
 
 (docu-units-scaling)=
 
@@ -90,7 +90,7 @@ The analytical solutions provide extreme performance. Single field evaluations t
 
 ## Numerical stability
 
-Many expressions provided in the literature have very questionable numerical stability. Most of these issues are fixed in Magpylib, but one should be aware that the results might not have more than 8-10 correct significant figures. A detailed treatsie on this topic is work in progress.
+Many expressions provided in the literature have very questionable numerical stability. Many of these problems are fixed in Magpylib, but one should be aware that accuracy can be a problem very close to objects, close the z-axis in cylindrical symmetries, at adge extensions, and at large distances. We are working on fixing these problems. Some details can be found in \[Ortner2022\].
 
 **References**
 
@@ -102,9 +102,8 @@ Many expressions provided in the literature have very questionable numerical sta
 - \[Derby2009\] N. Derby, "Cylindrical Magnets and Ideal Solenoids", arXiv:0909.3880v1, 2009
 - \[Caciagli2018\] A. Caciagli, R. J. Baars, A. P. Philipse and B. W. M. Kuipers, "Exact expression for the magnetic field of a finite cylinder with arbitrary uniform magnetization", Journal of Magnetism and Magnetic Materials 456 (2018) 423–432.
 - \[Slanovc2022\] F. Slanovc, M. Ortner, M. Moridi, C. Abert and D. Suess, "Full analytical solution for the magnetic field of uniformly magnetized cylinder tiles", submitted to Journal of Magnetism and Magnetic Materials.
-- \[Smythe1950\] W.B. Smythe, "Static and dynamic electricity" McGraw-Hill New York, 1950, vol. 3.
-- \[Simpson2001\] J. Simplson et al., "Simple analytic expressions for the magnetic field of a circular current loop," 2001.
-- \[Ortner2017\] M. Ortner et al., "Feedback of Eddy Currents in Layered Materials for Magnetic Speed Sensing", IEEE Transactions on Magnetics ( Volume: 53, Issue: 8, Aug. 2017)
+- \[Ortner2022\] M. Ortner, S. Slanovc and P. Leitner, "Numerically Stable and Computationally Efficient Expression for the Magnetic Field of a Current Loop", MDPI Magnetism, 3(1), 11-31, 2022.
+- \[Guptasarma1911\] D. Guptasarma and B. Singh, "New scheme for computing the magnetic field resulting from a uniformly magnetized arbitrary polyhedron", Geophysics (1999), 64(1):70.
 - \[Janssen2009\] J.L.G. Janssen, J.J.H. Paulides and E.A. Lomonova, "3D ANALYTICAL FIELD CALCULATION USING TRIANGULAR MAGNET SEGMENTS APPLIED TO A SKEWED LINEAR PERMANENT MAGNET ACTUATOR", ISEF 2009 - XIV International Symposium on Electromagnetic Fields in Mechatronics, Electrical and Electronic Engineering Arras, France, September 10-12, 2009
 - \[Rubeck2013\] C. Rubeck et al., "Analytical Calculation of Magnet Systems: Magnetic Field Created by Charged Triangles and Polyhedra", IEEE Transactions on Magnetics, VOL. 49, NO. 1, 2013
 - \[Jackson1999\] J. D. Jackson, "Classical Electrodynamics", 1999 Wiley, New York
diff --git a/docs/examples/01_fundamentals.md b/docs/examples/01_fundamentals.md
deleted file mode 100644
index be376f176..000000000
--- a/docs/examples/01_fundamentals.md
+++ /dev/null
@@ -1,28 +0,0 @@
----
-jupytext:
-  text_representation:
-    extension: .md
-    format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.13.7
-kernelspec:
-  display_name: Python 3 (ipykernel)
-  language: python
-  name: python3
----
-
-This section includes code examples. Detailed package, class, method and function documentations are found in the {ref}`modindex`. A technical introduction is given in {ref}`intro`.
-
-Fundamentals
-============
-
-```{toctree}
-:caption: Fundamentals
-:maxdepth: 2
-
-examples_02_paths.md
-examples_05_backend_canvas.md
-examples_03_collections.md
-examples_01_complex_forms.md
-examples_04_custom_source.md
-```
diff --git a/docs/examples/02_graphic_output.md b/docs/examples/02_graphic_output.md
deleted file mode 100644
index 949b6f185..000000000
--- a/docs/examples/02_graphic_output.md
+++ /dev/null
@@ -1,29 +0,0 @@
----
-jupytext:
-  text_representation:
-    extension: .md
-    format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.13.1
-kernelspec:
-  display_name: Python 3
-  language: python
-  name: python3
----
-
-
-This section includes code examples. Detailed package, class, method and function documentations are found in the {ref}`modindex`. A technical introduction is given in {ref}`intro`.
-
-(example-gallery-graphic-output)=
-
-Graphic Output
-==============
-
-```{toctree}
-:caption: Graphic Output
-:maxdepth: 2
-
-examples_12_styles.md
-examples_13_3d_models.md
-examples_15_animation.md
-```
\ No newline at end of file
diff --git a/docs/examples/03_advanced.md b/docs/examples/03_advanced.md
deleted file mode 100644
index 41f1ba1f1..000000000
--- a/docs/examples/03_advanced.md
+++ /dev/null
@@ -1,25 +0,0 @@
----
-jupytext:
-  text_representation:
-    extension: .md
-    format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.13.1
-kernelspec:
-  display_name: Python 3
-  language: python
-  name: python3
----
-
-This section includes code examples. Detailed package, class, method and function documentations are found in the {ref}`modindex`. A technical introduction is given in {ref}`intro`.
-
-Advanced
-========
-
-```{toctree}
-:caption: Advanced Features
-:maxdepth: 2
-
-examples_21_compound.md
-examples_22_field_interpolation.md
-```
\ No newline at end of file
diff --git a/docs/examples/04_application_examples.md b/docs/examples/04_application_examples.md
deleted file mode 100644
index 48312bee2..000000000
--- a/docs/examples/04_application_examples.md
+++ /dev/null
@@ -1,25 +0,0 @@
----
-jupytext:
-  text_representation:
-    extension: .md
-    format_name: myst
-    format_version: 0.13
-    jupytext_version: 1.13.1
-kernelspec:
-  display_name: Python 3
-  language: python
-  name: python3
----
-
-This section includes code examples. Detailed package, class, method and function documentations are found in the {ref}`modindex`. A technical introduction is given in {ref}`intro`.
-
-Application Examples
-====================
-
-```{toctree}
-:caption: Application Examples
-:maxdepth: 2
-
-examples_30_coil_field_lines.md
-examples_31_end_of_shaft.md
-```
\ No newline at end of file
diff --git a/docs/examples/examples_01_complex_forms.md b/docs/examples/examples_01_complex_forms.md
index dfb822102..411fabef6 100644
--- a/docs/examples/examples_01_complex_forms.md
+++ b/docs/examples/examples_01_complex_forms.md
@@ -13,9 +13,9 @@ kernelspec:
 
 (examples-complex-forms)=
 
-# Complex forms
+# Complex shapes - Superposition
 
-The [**superposition principle**](https://en.wikipedia.org/wiki/Superposition_principle) states that the net response caused by two or more stimuli is the sum of the responses caused by each stimulus individually. This principle holds in Magneto statics when there is no material response, and simply means that the total field created by multiple magnets is the sum of the individual fields.
+The [**superposition principle**](https://en.wikipedia.org/wiki/Superposition_principle) states that the net response caused by two or more stimuli is the sum of the responses caused by each stimulus individually. This principle holds in Magneto statics when there is no material response, and simply means that the total field created by multiple magnets and currents is the sum of the individual fields.
 
 It is critical to understand that the superposition principle holds for the magnetization itself. When two magnets overlap geometrically, the magnetization in the overlap region is given by the vector sum of the two individual magnetizations.
 
@@ -84,4 +84,6 @@ print('getB from Cylindersegment', ring0.getB((1,2,3)))
 print('getB from Cylinder cut-out', ring1.getB((1,2,3)))
 ```
 
-Note that, it is faster to compute the `Cylinder` field two times than computing the complex `CylinderSegment` field one time. This is why Magpylib automatically falls back to the `Cylinder` solution whenever `CylinderSegment` is called with 360 deg section angles. Unfortunately, cut-out operations cannot be displayed graphically at the moment, but {ref}`examples-own-3d-models` offer a solution here.
\ No newline at end of file
+Note that, it is faster to compute the `Cylinder` field two times than computing the complex `CylinderSegment` field one time. This is why Magpylib automatically falls back to the `Cylinder` solution whenever `CylinderSegment` is called with 360 deg section angles. Unfortunately, cut-out operations cannot be displayed graphically at the moment, but {ref}`examples-own-3d-models` offer a solution here.
+
+Finally, it is explained in {ref}`examples-triangle`, how complex shapes are achieved based on triangular meshes.
diff --git a/docs/examples/examples_02_paths.md b/docs/examples/examples_02_paths.md
index b448ea7bd..366620b60 100644
--- a/docs/examples/examples_02_paths.md
+++ b/docs/examples/examples_02_paths.md
@@ -15,6 +15,9 @@ kernelspec:
 
 # Paths
 
+Magpylib objects all have `position` and `orientation` attributes in the global reference frame. These attributes can be single instances describing a single position and a single orientation of an object, but they can also be vectors of positions and orientations. In this case the attributes represent multiple object locations and orientaions and are referred to as a **path**. When the field of an object is computed, it is automatically computed for the whole path. The reason is that computation is much more efficient when all instances are computed in a single operation. The introduction to `position` and `orientation` attributes and path can be found in {ref}`intro-position-and-orientation`. Here we show how to work with paths properly.
+
+
 (examples-assign-absolute-path)=
 ## Assigning absolute paths
 
diff --git a/docs/examples/examples_03_collections.md b/docs/examples/examples_03_collections.md
index 984059669..82b471bed 100644
--- a/docs/examples/examples_03_collections.md
+++ b/docs/examples/examples_03_collections.md
@@ -15,9 +15,11 @@ kernelspec:
 
 # Collections
 
+The `Collection` class is a versatile way of grouping and manipulating multiple Magpylib objects. A basic introduction is given in {ref}`intro-collections`. Here things are explained in more detail with examples.
+
 ## Constructing collections
 
-The `Collection` class is a versatile way of grouping and manipulating Magpylib objects. When objects are added to a Collection they are added by reference (not copied) to the **attributes** `children` (list of all objects), `sources` (list of the sources), `sensors` (list of the sensors) and `collections`.
+Collections have the attributes `children`, `sources`, `sensors` and `collections`. These attributes are ordered lists that contain objects that are added to the collection by reference (not copied). `children` returns is list of all objects in the collection. `sources` returns a list of the sources, `sensors` a list of the sensors and `collections` a list of "sub-collections" within the collection.
 
 ```{code-cell} ipython3
 import magpylib as magpy
@@ -34,7 +36,7 @@ print(f"sensors:     {coll.sensors}")
 print(f"collections: {coll.collections}")
 ```
 
-These attributes are ordered lists. New additions are always added at the end. Add objects to an existing collection using these parameters, or the **`add`** method.
+New additions are always added at the end. Add objects to an existing collection using these parameters, or the **`add`** method.
 
 ```{code-cell} ipython3
 
diff --git a/docs/examples/examples_04_custom_source.md b/docs/examples/examples_04_custom_source.md
index 613a6889a..aa758b9f3 100644
--- a/docs/examples/examples_04_custom_source.md
+++ b/docs/examples/examples_04_custom_source.md
@@ -12,7 +12,7 @@ kernelspec:
 ---
 
 (examples-custom-source-objects)=
-# Custom source objects
+# CustomSource class
 
 The class `CustomSource` allows users to integrate their own custom-object field computation into the Magpylib interface. For this, the argument `field_func` must be provided with a function that is then automatically called with `getB` and `getH`. This custom field function is treated like a core function. It must have the positional arguments `field` (with values `'B'` or `'H'`), and `observers` (must accept array_like, shape (n,3)) and return the B-field in units of \[mT\] and the H-field in units of \[kA/m\] with a similar shape. A fundamental example how to create a custom source object is shown below:
 
diff --git a/docs/examples/examples_05_backend_canvas.md b/docs/examples/examples_05_backend_canvas.md
index d12c784d8..ddf18705a 100644
--- a/docs/examples/examples_05_backend_canvas.md
+++ b/docs/examples/examples_05_backend_canvas.md
@@ -13,7 +13,9 @@ kernelspec:
 
 (examples-backends-canvas)=
 
-# Backend and canvas
+# Graphics - Backend and canvas
+
+The graphic backend refers to the plotting library that is used for graphic output. Canvas refers to the frame/window/canvas/axes object the graphic output is forwarded to.
 
 ## Graphic backend
 
diff --git a/docs/examples/examples_06_triangle.md b/docs/examples/examples_06_triangle.md
new file mode 100644
index 000000000..1db4dfadb
--- /dev/null
+++ b/docs/examples/examples_06_triangle.md
@@ -0,0 +1,106 @@
+---
+jupytext:
+  text_representation:
+    extension: .md
+    format_name: myst
+    format_version: 0.13
+    jupytext_version: 1.13.7
+kernelspec:
+  display_name: Python 3 (ipykernel)
+  language: python
+  name: python3
+---
+
+(examples-triangle)=
+
+# Complex shapes - The Triangle class
+
+The field of a homogeneously charged magnet is, on the outside, the same as the field of a similarly shaped body with a magnetic surface charge. The surface charge is proportional to the projection of the magnetization vector onto the surface normal. The `Triangle` class is set up so that it can easily be used to approximate surfaces, and given the magnetization vector, the charge density is automatically computed. The resulting H-field is correct, but the B-field is only correct on the outside of the body, because on the inside the magnetization must be added to it.
+
+## Example: Triangular prisma magnet
+
+Consider a Prisma with triangular base that is magnetized orthogonal to the base. All surface normals of the sides of the prisma are orthogonal to the magnetization vector. As a result the sides do not contribute to the magnetic field because their charge desity dissappears. Only top and bottom surfaces contribute. One must be very careful when defining those surfaces in such a way that the surface normals point outwards. The following examples shows how the `Triangle` class can be used to create such a prisma.
+
+```{code-cell} ipython3
+import magpylib as magpy
+
+top = magpy.misc.Triangle(
+    magnetization=(0,0,1000),
+    vertices= ((-1,-1,1), (1,-1,1), (0,2,1)),
+)
+bott = magpy.misc.Triangle(
+    magnetization=(0,0,1000),
+    vertices= ((-1,-1,-1), (0,2,-1), (1,-1,-1)),
+)
+
+prisma = magpy.Collection(top, bott)
+
+prisma.show(style_magnetization_size=0.2)
+```
+
+## Example: Cuboctahedron magnet
+
+More complex bodies are easy constructed from Triangles. The following code shows how a magnet with cuboctahedron shape can be constructed. Be aware that the B-field is only correcto on the outside.
+
+```{code-cell} ipython3
+import magpylib as magpy
+
+vertices = (((0.0,1.0,-1.0),(1.0,1.0,0.0),(-1.0,1.0,0.0)),
+    ((0.0,1.0,1.0),(-1.0,1.0,0.0),(1.0,1.0,0.0)),
+    ((0.0,1.0,1.0),(0.0,-1.0,1.0),(-1.0,0.0,1.0)),
+    ((0.0,1.0,1.0),(1.0,0.0,1.0),(0.0,-1.0,1.0)),
+    ((0.0,1.0,-1.0),(0.0,-1.0,-1.0),(1.0,0.0,-1.0)),
+    ((0.0,1.0,-1.0),(-1.0,0.0,-1.0),(0.0,-1.0,-1.0)),
+    ((0.0,-1.0,1.0),(1.0,-1.0,0.0),(-1.0,-1.0,0.0)),
+    ((0.0,-1.0,-1.0),(-1.0,-1.0,0.0),(1.0,-1.0,0.0)),
+    ((-1.0,1.0,0.0),(-1.0,0.0,1.0),(-1.0,0.0,-1.0)),
+    ((-1.0,-1.0,0.0),(-1.0,0.0,-1.0),(-1.0,0.0,1.0)),
+    ((1.0,1.0,0.0),(1.0,0.0,-1.0),(1.0,0.0,1.0)),
+    ((1.0,-1.0,0.0),(1.0,0.0,1.0),(1.0,0.0,-1.0)),
+    ((0.0,1.0,1.0),(-1.0,0.0,1.0),(-1.0,1.0,0.0)),
+    ((0.0,1.0,1.0),(1.0,1.0,0.0),(1.0,0.0,1.0)),
+    ((0.0,1.0,-1.0),(-1.0,1.0,0.0),(-1.0,0.0,-1.0)),
+    ((0.0,1.0,-1.0),(1.0,0.0,-1.0),(1.0,1.0,0.0)),
+    ((0.0,-1.0,-1.0),(-1.0,0.0,-1.0),(-1.0,-1.0,0.0)),
+    ((0.0,-1.0,-1.0),(1.0,-1.0,0.0),(1.0,0.0,-1.0)),
+    ((0.0,-1.0,1.0),(-1.0,-1.0,0.0),(-1.0,0.0,1.0)),
+    ((0.0,-1.0,1.0),(1.0,0.0,1.0),(1.0,-1.0,0.0)),
+)
+cubo = magpy.Collection()
+for vert in vertices:
+    cubo.add(magpy.misc.Triangle(
+        magnetization=(100,200,300),
+        vertices=vert,
+    ))
+
+magpy.show(*cubo)
+```
+
+
+<!-- ## Pyvista mesh and Facet class # we leave this for the TriMesh class :)
+
+Contemporary tools like [Pyvista](https://docs.pyvista.org/) offer powerful meshing options. The following example shows how a complex Pyvista object can be used together with the `Facet` class to create a magnet object with little effort. However, contrary to a an actual magnet object there is no insode-outside check when transofming between B- and H- field.
+
+```{code-cell} ipython3
+import pyvista as pv
+import magpylib as magpy
+
+# create a complex pyvista PolyData object
+sphere = pv.Sphere(radius=0.85)
+dodec = pv.Dodecahedron().triangulate().subdivide(5)
+object = dodec.boolean_difference(sphere)
+
+# extract triangles and create Facet vertices input
+points = object.points
+faces = object.faces.reshape(-1,4)[:,1:]
+verts = [[points[f[i]] for i in range(3)] for f in faces]
+
+magnet = magpy.misc.Facet(
+    magnetization=(0,0,100),
+    vertices=verts
+)
+
+magnet.show(backend='plotly')
+``` -->
+
+
diff --git a/docs/examples/examples_12_styles.md b/docs/examples/examples_12_styles.md
index 9fb8fa407..78e0eeffc 100644
--- a/docs/examples/examples_12_styles.md
+++ b/docs/examples/examples_12_styles.md
@@ -12,7 +12,7 @@ kernelspec:
 ---
 
 (examples-graphic-styles)=
-# Styles
+# Graphics - Styles
 
 The graphic styles define how Magpylib objects are displayed visually when calling `show`. They can be fine-tuned and individualized in many ways.
 
diff --git a/docs/examples/examples_13_3d_models.md b/docs/examples/examples_13_3d_models.md
index 057c0dd07..dcced61e8 100644
--- a/docs/examples/examples_13_3d_models.md
+++ b/docs/examples/examples_13_3d_models.md
@@ -13,7 +13,7 @@ kernelspec:
 
 (examples-3d-models)=
 
-# 3D models and CAD
+# Graphics - 3D models and CAD
 
 (examples-own-3d-models)=
 ## Custom 3D models
diff --git a/docs/examples/examples_15_animation.md b/docs/examples/examples_15_animation.md
index c2b052684..167f32ec2 100644
--- a/docs/examples/examples_15_animation.md
+++ b/docs/examples/examples_15_animation.md
@@ -13,9 +13,9 @@ kernelspec:
 
 (examples-animation)=
 
-# Animate paths
+# Graphics - Animate paths
 
-With the Plotly backend, paths can automatically be animated with `show(animation=True)`. Animations can be fine-tuned with the following properties:
+With some backends, paths can automatically be animated with `show(animation=True)`. Animations can be fine-tuned with the following properties:
 
 1. `animation_time` (default=3), must be a positive number that gives the animation time in seconds.
 2. `animation_slider` (default=`True`), is boolean and sets if a slider should be displayed in addition.
diff --git a/docs/examples/examples_21_compound.md b/docs/examples/examples_21_compound.md
index 4f1053f72..6b9369dcc 100644
--- a/docs/examples/examples_21_compound.md
+++ b/docs/examples/examples_21_compound.md
@@ -13,7 +13,7 @@ kernelspec:
 
 (examples-compounds)=
 
-# Compounds
+# Collections - Compounds
 
 The `Collection` class is a powerful tool for grouping and tracking object assemblies.
 However, in many cases it is convenient to have assembly variables themselves (e.g. geometric arrangement) as class properties of new custom classes, which is achieved by sub-classing `Collection`. We refer to such super-classes as **compounds** and show how to seamlessly integrate them into Magpylib.
diff --git a/docs/examples/examples_30_coil_field_lines.md b/docs/examples/examples_30_coil_field_lines.md
index a799e1a58..48e048b91 100644
--- a/docs/examples/examples_30_coil_field_lines.md
+++ b/docs/examples/examples_30_coil_field_lines.md
@@ -11,7 +11,7 @@ kernelspec:
   name: python3
 ---
 
-# Coil Field Lines
+# Magnetic field lines
 
 In this example we model the **magnetic field of a coil**, and show how to display it with spectacular **field line** representations.
 
diff --git a/docs/examples/examples_31_end_of_shaft.md b/docs/examples/examples_31_end_of_shaft.md
index 2988514d2..1609f68ca 100644
--- a/docs/examples/examples_31_end_of_shaft.md
+++ b/docs/examples/examples_31_end_of_shaft.md
@@ -11,7 +11,7 @@ kernelspec:
   name: python3
 ---
 
-# End of shaft angle sensor
+# Angle sensor
 
 End of shaft angle sensing is a classical example for a magnetic position system. The goal is to determine the angular position of a rotating shaft. A magnet, typically a diametrically magnetized cylinder, is mounted at the end of the shaft. A 2D sensor is mounted below. When the shaft rotates the two sensor outputs will be $s_1=B_0 sin(\varphi)$ and $s_2=B_0 cos(\varphi)$, so that the angle is uniquely given by $\varphi = arctan(s_1/s_2)$.
 
@@ -19,27 +19,26 @@ In the example below we show such a typical end-of-shaft system with a 2-pixel s
 
 ```{code-cell} ipython3
 import numpy as np
-import plotly.graph_objects as go
+import plotly.express as px
 import magpylib as magpy
 
 # create magnet
 magnet = magpy.magnet.Cylinder(
-    magnetization=(1000,0,0),
-    dimension=(6,2),
-    position=(0,0,1.5),
-    style_label='Magnet',
-    style_color='.7',
+    magnetization=(1000, 0, 0),
+    dimension=(6, 2),
+    position=(0, 0, 1.5),
+    style_label="Magnet",
+    style_color=".7",
 )
 
 # create shaft dummy with 3D model
 shaft = magpy.misc.CustomSource(
-    position=(0,0,7),
-    style_color='.7',
+    position=(0, 0, 7),
+    style_color=".7",
     style_model3d_showdefault=False,
-    style_label='Shaft'
+    style_label="Shaft",
 )
 shaft_trace = magpy.graphics.model3d.make_Prism(
-    backend='plotly',
     base=20,
     diameter=10,
     height=10,
@@ -52,35 +51,34 @@ displacement = 1
 angles = np.linspace(0, 360, 72)
 coll = magnet + shaft
 magnet.move((displacement, 0, 0))
-coll.rotate_from_angax(angles, 'z', anchor=0, start=0)
+coll.rotate_from_angax(angles, "z", anchor=0, start=0)
 
 # create sensor
 gap = 3
 sens = magpy.Sensor(
-    position=(0,0,-gap),
-    pixel=[(1,0,0), (-1,0,0)],
+    position=(0, 0, -gap),
+    pixel=[(1, 0, 0), (-1, 0, 0)],
     style_pixel_size=0.5,
     style_size=1.5,
 )
 
 # show 3D animation of wobble motion
 fig1 = go.Figure()
-magpy.show(magnet, sens, shaft, animation=True, backend='plotly', canvas=fig1)
+magpy.show(magnet, sens, shaft, animation=True, backend="plotly", canvas=fig1)
 fig1.update_layout(scene_camera_eye_z=-1.1)
 fig1.show()
 
 # show sensor output in plotly
 fig2 = go.Figure()
-B = sens.getB(magnet)
-for px,dash in zip([0,1], ['solid', 'dash']):
-    for i,xy,col in zip([0,1], ['x','y'], ['red', 'green']):
-        fig2.add_trace(go.Scatter(x=angles, y=B[:,px,i],
-            name=f"pixel{px}-B{xy}",
-            line=dict(color=col, dash=dash),
-        ))
-fig2.update_layout(
-    xaxis_title='angle [deg]',
-    yaxis_title='field [mT]'
+df = sens.getB(magnet, output="dataframe")
+df["angle [deg]"] = angles[df["path"]]
+
+fig2 = px.line(
+    df,
+    x="angle [deg]",
+    y=["Bx", "By"],
+    line_dash="pixel",
+    labels={"value": "Field [mT]"},
 )
 fig2.show()
 ```
diff --git a/docs/index.md b/docs/index.md
index bf5ddff7f..73251c08c 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -5,7 +5,7 @@
 
 # What is Magpylib ?
 
-Magpylib is a Python package for calculating **3D static magnetic fields** of magnets, line currents and other sources. The computation is based on analytical expressions and therefore **extremely fast**. A **user friendly geometry interface** enables convenient relative positioning between sources and observers.
+Magpylib is a Python package for calculating **3D static magnetic fields** of magnets, line currents and other sources. The computation is based on analytical expressions and therefore **extremely fast**. A **user friendly geometry API** enables convenient relative positioning between sources and observers.
 
 ```{image} _static/images/index/source_fundamentals.png
 :alt: source fundamentals
@@ -42,15 +42,31 @@ Here, a cylinder shaped permanent magnet with (diameter, height) of (4, 5) milli
 _pages/*
 ```
 
+```{toctree}
+:caption: FUNCTIONALITY
+:glob: true
+:maxdepth: 2
+
+examples/examples_02_paths.md
+examples/examples_05_backend_canvas.md
+examples/examples_12_styles.md
+examples/examples_15_animation.md
+examples/examples_13_3d_models.md
+examples/examples_03_collections.md
+examples/examples_21_compound.md
+examples/examples_01_complex_forms.md
+examples/examples_06_triangle.md
+examples/examples_04_custom_source.md
+```
+
 ```{toctree}
 :caption: EXAMPLE GALLERIES
 :glob: true
 :maxdepth: 2
 
-examples/01_fundamentals.md
-examples/02_graphic_output.md
-examples/03_advanced.md
-examples/04_application_examples.md
+examples/examples_30_coil_field_lines.md
+examples/examples_22_field_interpolation.md
+examples/examples_31_end_of_shaft.md
 ```
 
 ```{toctree}
diff --git a/magpylib/__init__.py b/magpylib/__init__.py
index 073676593..93aafdc21 100644
--- a/magpylib/__init__.py
+++ b/magpylib/__init__.py
@@ -28,7 +28,7 @@
 
 """
 # module level dunders
-__version__ = "4.1.2"
+__version__ = "4.2.0"
 __author__ = "Michael Ortner & Alexandre Boisselet"
 __credits__ = "The Magpylib community"
 __all__ = [
diff --git a/magpylib/_src/display/traces_base.py b/magpylib/_src/display/traces_base.py
index 44fed8bc9..770a9ecf6 100644
--- a/magpylib/_src/display/traces_base.py
+++ b/magpylib/_src/display/traces_base.py
@@ -2,6 +2,7 @@
 from functools import partial
 
 import numpy as np
+from scipy.spatial import ConvexHull  # pylint: disable=no-name-in-module
 
 from magpylib._src.display.traces_utility import merge_mesh3d
 from magpylib._src.display.traces_utility import place_and_orient_model3d
@@ -532,3 +533,138 @@ def make_Arrow(
     trace = merge_mesh3d(cone, prism)
     trace = place_and_orient_model3d(trace, orientation=orientation, position=position)
     return get_model(trace, backend=backend, show=show, scale=scale, kwargs=kwargs)
+
+
+def make_Tetrahedron(
+    backend="generic",
+    vertices=None,
+    position=None,
+    orientation=None,
+    show=True,
+    scale=1,
+    **kwargs,
+) -> dict:
+    """Provides the 3D-model parameters for a pyramid in dictionary form, based on
+    number of vertices of the base, diameter and height. The zero position is in the
+    barycenter of the vertices.
+
+    Parameters
+    ----------
+    backend : str
+        Plotting backend corresponding to the trace. Can be one of `['matplotlib', 'plotly']`.
+
+    vertices: ndarray, shape (4,3)
+        Vertices (x1,y1,z1), (x2,y2,z2), (x3,y3,z3), (x4,y4,z4), in the relative
+        coordinate system of the tetrahedron.
+
+    position : array_like, shape (3,), default=(0,0,0)
+        Reference position of the vertices in the global CS. The zero position is
+        in the barycenter of the vertices.
+
+    orientation : scipy Rotation object with length 1 or m, default=`identity`
+        Orientation of the vertices in the global CS.
+
+    show : bool, default=True
+        Shows/hides model3d object based on provided trace.
+
+    scale : float, default=1
+        Scaling factor by which the trace vertices coordinates are multiplied.
+
+    **kwargs : dict, optional
+        Additional keyword arguments to be passed to the trace constructor directly.
+        (e.g. `opacity=0.5` for plotly or `alpha=0.5` for matplotlib)
+
+    Returns
+    -------
+    3D-model: dict
+        A dictionary with necessary key/value pairs with the necessary information to construct
+        a 3D-model.
+    """
+    x, y, z = np.array(vertices).T
+    trace = dict(
+        i=np.array(
+            [
+                0,
+                0,
+                1,
+                2,
+            ]
+        ),
+        j=np.array([1, 1, 2, 0]),
+        k=np.array([2, 3, 3, 3]),
+        x=x,
+        y=y,
+        z=z,
+    )
+
+    trace = place_and_orient_model3d(trace, orientation=orientation, position=position)
+    return get_model(trace, backend=backend, show=show, scale=scale, kwargs=kwargs)
+
+
+def make_TriangularMesh(
+    backend="generic",
+    vertices=None,
+    triangles=None,
+    position=None,
+    orientation=None,
+    show=True,
+    scale=1,
+    **kwargs,
+) -> dict:
+    """Provides the 3D-model parameters for a custom triangular mesh in dictionary form, based on
+    number of vertices of the base, diameter and height. The zero position is in the
+    barycenter of the vertices.
+
+    Parameters
+    ----------
+    backend : str
+        Plotting backend corresponding to the trace. Can be one of `['matplotlib', 'plotly']`.
+
+    vertices: ndarray, shape (4,3)
+        Vertices (x1,y1,z1), (x2,y2,z2), (x3,y3,z3), (x4,y4,z4), in the relative
+        coordinate system of the triangular mesh.
+
+    triangles: ndarray, shape (4,3)
+        For each triangle, the indices of the three points that make up the triangle, ordered in an
+        anticlockwise manner. If not specified, a `scipy.spatial.ConvexHull` triangulation is
+        calculated.
+
+    position : array_like, shape (3,), default=(0,0,0)
+        Reference position of the vertices in the global CS. The zero position is
+        in the barycenter of the vertices.
+
+    orientation : scipy Rotation object with length 1 or m, default=`identity`
+        Orientation of the vertices in the global CS.
+
+    show : bool, default=True
+        Shows/hides model3d object based on provided trace.
+
+    scale : float, default=1
+        Scaling factor by which the trace vertices coordinates are multiplied.
+
+    **kwargs : dict, optional
+        Additional keyword arguments to be passed to the trace constructor directly.
+        (e.g. `opacity=0.5` for plotly or `alpha=0.5` for matplotlib)
+
+    Returns
+    -------
+    3D-model: dict
+        A dictionary with necessary key/value pairs with the necessary information to construct
+        a 3D-model.
+    """
+    vertices = np.array(vertices)
+    x, y, z = vertices.T
+    if triangles is None:
+        hull = ConvexHull(vertices)
+        triangles = hull.simplices
+    i, j, k = np.array(triangles).T
+    trace = dict(
+        i=i,
+        j=j,
+        k=k,
+        x=x,
+        y=y,
+        z=z,
+    )
+    trace = place_and_orient_model3d(trace, orientation=orientation, position=position)
+    return get_model(trace, backend=backend, show=show, scale=scale, kwargs=kwargs)
diff --git a/magpylib/_src/display/traces_generic.py b/magpylib/_src/display/traces_generic.py
index faf7e1b59..09a94f7fe 100644
--- a/magpylib/_src/display/traces_generic.py
+++ b/magpylib/_src/display/traces_generic.py
@@ -10,6 +10,7 @@
 import numpy as np
 from scipy.spatial.transform import Rotation as RotScipy
 
+import magpylib as magpy
 from magpylib import _src
 from magpylib._src.defaults.defaults_classes import default_settings as Config
 from magpylib._src.defaults.defaults_utility import linearize_dict
@@ -21,6 +22,10 @@
 )
 from magpylib._src.display.traces_base import make_Ellipsoid as make_BaseEllipsoid
 from magpylib._src.display.traces_base import make_Prism as make_BasePrism
+from magpylib._src.display.traces_base import make_Tetrahedron as make_BaseTetrahedron
+from magpylib._src.display.traces_base import (
+    make_TriangularMesh as make_BaseTriangularMesh,
+)
 from magpylib._src.display.traces_utility import draw_arrow_from_vertices
 from magpylib._src.display.traces_utility import draw_arrowed_circle
 from magpylib._src.display.traces_utility import draw_arrowed_line
@@ -41,7 +46,7 @@
 
 
 class MagpyMarkers:
-    """A class that stores markers 3D-coordinates"""
+    """A class that stores markers 3D-coordinates."""
 
     def __init__(self, *markers):
         self.style = Markers()
@@ -345,6 +350,60 @@ def make_Sphere(
     )
 
 
+def make_Tetrahedron(
+    obj,
+    position=(0.0, 0.0, 0.0),
+    orientation=None,
+    color=None,
+    style=None,
+    **kwargs,
+) -> dict:
+    """
+    Create the plotly mesh3d parameters for a Tetrahedron Magnet in a dictionary based on the
+    provided arguments.
+    """
+    style = obj.style if style is None else style
+    trace = make_BaseTetrahedron("plotly-dict", vertices=obj.vertices, color=color)
+    update_trace_name(trace, "Tetrahedron", "", style)
+    update_magnet_mesh(
+        trace, mag_style=style.magnetization, magnetization=obj.magnetization
+    )
+    return place_and_orient_model3d(
+        trace, orientation=orientation, position=position, **kwargs
+    )
+
+
+def make_Triangle(
+    obj,
+    position=(0.0, 0.0, 0.0),
+    orientation=None,
+    color=None,
+    style=None,
+    **kwargs,
+) -> dict:
+    """
+    Creates the plotly mesh3d parameters for a TriangularMesh Magnet in a dictionary based on the
+    provided arguments.
+    """
+    vertices = obj.vertices
+
+    # Return unique vertices and triangles from vertices
+    vert_unique, tr = np.unique(vertices.reshape((-1, 3)), axis=0, return_inverse=True)
+    triangles = tr.reshape((-1, 3))
+
+    style = obj.style if style is None else style
+    trace = make_BaseTriangularMesh(
+        "plotly-dict", vertices=vert_unique, triangles=triangles, color=color
+    )
+    update_trace_name(trace, "Triangle", "", style)
+    update_magnet_mesh(
+        trace, mag_style=style.magnetization, magnetization=obj.magnetization
+    )
+    return place_and_orient_model3d(
+        trace, orientation=orientation, position=position, **kwargs
+    )
+
+
 def make_Pixels(positions, size=1) -> dict:
     """
     Create the plotly mesh3d parameters for Sensor pixels based on pixel positions and chosen size
@@ -490,6 +549,10 @@ def make_mag_arrows(obj, style, legendgroup, kwargs):
     # vector length, color and magnetization
     if hasattr(obj, "diameter"):
         length = obj.diameter  # Sphere
+    elif isinstance(obj, magpy.misc.Triangle):
+        length = np.amax(obj.vertices) - np.amin(obj.vertices)
+    elif hasattr(obj, "vertices"):
+        length = np.amax(np.ptp(obj.vertices, axis=0))
     else:  # Cuboid, Cylinder, CylinderSegment
         length = np.amax(obj.dimension[:3])
     length *= 1.8 * style.magnetization.size
diff --git a/magpylib/_src/fields/field_BH_cylinder.py b/magpylib/_src/fields/field_BH_cylinder.py
index 721cf957a..d695c47e8 100644
--- a/magpylib/_src/fields/field_BH_cylinder.py
+++ b/magpylib/_src/fields/field_BH_cylinder.py
@@ -247,7 +247,7 @@ def fieldH_cylinder_diametral(
     return Hr, Hphi, Hz
 
 
-# ON INTERFACE
+# CORE
 def magnet_cylinder_field(
     field: str,
     observers: np.ndarray,
diff --git a/magpylib/_src/fields/field_BH_cylinder_segment.py b/magpylib/_src/fields/field_BH_cylinder_segment.py
index 4f49b39bf..466149a29 100644
--- a/magpylib/_src/fields/field_BH_cylinder_segment.py
+++ b/magpylib/_src/fields/field_BH_cylinder_segment.py
@@ -2345,6 +2345,7 @@ def magnet_cylinder_segment_field_internal(
     return BHfinal
 
 
+# CORE
 def magnet_cylinder_segment_field(
     field: str,
     observers: np.ndarray,
diff --git a/magpylib/_src/fields/field_BH_dipole.py b/magpylib/_src/fields/field_BH_dipole.py
index 6b6e25f85..db259359b 100644
--- a/magpylib/_src/fields/field_BH_dipole.py
+++ b/magpylib/_src/fields/field_BH_dipole.py
@@ -6,6 +6,7 @@
 from magpylib._src.input_checks import check_field_input
 
 
+# CORE
 def dipole_field(
     field: str,
     observers: np.ndarray,
@@ -19,10 +20,10 @@ def dipole_field(
     ----------
     field: str, default=`'B'`
         If `field='B'` return B-field in units of [mT], if `field='H'` return H-field
+        in units of [kA/m].
 
     observers: ndarray, shape (n,3)
         Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
-        in units of [kA/m].
 
     moment: ndarray, shape (n,3)
         Dipole moment vector in units of [mT*mm^3].
@@ -34,6 +35,8 @@ def dipole_field(
 
     Examples
     --------
+    Compute the B-field of two different dipole-observer instances.
+
     >>> import numpy as np
     >>> import magpylib as magpy
     >>> mom = np.array([(1,2,3), (0,0,1)])
diff --git a/magpylib/_src/fields/field_BH_line.py b/magpylib/_src/fields/field_BH_line.py
index 92c8fa72e..10add772b 100644
--- a/magpylib/_src/fields/field_BH_line.py
+++ b/magpylib/_src/fields/field_BH_line.py
@@ -58,7 +58,7 @@ def current_vertices_field(
     return BH
 
 
-# ON INTERFACE
+# CORE
 def current_line_field(
     field: str,
     observers: np.ndarray,
diff --git a/magpylib/_src/fields/field_BH_loop.py b/magpylib/_src/fields/field_BH_loop.py
index 65a483007..cc76e5fc2 100644
--- a/magpylib/_src/fields/field_BH_loop.py
+++ b/magpylib/_src/fields/field_BH_loop.py
@@ -2,16 +2,15 @@
 Implementations of analytical expressions for the magnetic field of
 a circular current loop. Computation details in function docstrings.
 """
-# pylint: disable=no-name-in-module # reason scipy.special.ellipe import
 import numpy as np
-from scipy.special import ellipe
 
-from magpylib._src.fields.special_cel import cel_loop_stable
+from magpylib._src.fields.special_cel import cel_iter
 from magpylib._src.input_checks import check_field_input
 from magpylib._src.utility import cart_to_cyl_coordinates
 from magpylib._src.utility import cyl_field_to_cart
 
-# ON INTERFACE
+
+# CORE
 def current_loop_field(
     field: str,
     observers: np.ndarray,
@@ -36,6 +35,7 @@ def current_loop_field(
 
     diameter: ndarray, shape (n,)
         Diameter of loop in units of [mm].
+
     Returns
     -------
     B-field or H-field: ndarray, shape (n,3)
@@ -58,72 +58,66 @@ def current_loop_field(
 
     Notes
     -----
-    This field can be obtained by direct application of the Biot-Savardt law.
-    Several sources in the literature provides these formulas:
-
-    Smythe, "Static and dynamic electricity" McGraw-Hill New York, 1950, vol. 3.
-
-    Simpson, "Simple analytic expressions for the magnetic field of a circular current
-    loop," 2001.
-
-    Ortner, "Feedback of Eddy Currents in Layered Materials for Magnetic Speed Sensing",
-    IEEE Transactions on Magnetics ( Volume: 53, Issue: 8, Aug. 2017).
-
-    New numerically stable implementation based on [Ortner/Leitner wip]. Based
-    thereon the field is computed with >12 digits precision everywhere. On the
-    loop the result is set to (0,0,0).
+    Implementation based on "Numerically stable and computationally efficinet expression for
+    the magnetic field of a current loop.", M.Ortner et al, Submitted to MDPI Magnetism, 2022
     """
 
     bh = check_field_input(field, "current_loop_field()")
 
     r, phi, z = cart_to_cyl_coordinates(observers)
-    rad = np.abs(diameter / 2)
-    n = len(rad)
-
-    # define masks for special cases
-    mask_radius0 = rad == 0
-    # rel pos deviation by 1e-15 to account for num errors (e.g. when rotating)
-    mask_on_loop = np.logical_and(abs(r - rad) < 1e-15 * rad, z == 0)
-    mask_general = ~np.logical_or(mask_radius0, mask_on_loop)
-
-    # collect general case inputs
-    r = r[mask_general]
-    z = z[mask_general]
-    rad = rad[mask_general]
-    nX = len(rad)
-
-    # express through ratios (make dimensionless, avoid large/small input values)
-    rb = r / rad
-    zb = z / rad
-
-    # pre-compute small quantities that might not be cached
-    z2 = zb**2
-    brack = z2 + (rb + 1) ** 2
-    k2 = 4 * rb / brack
-    xi = cel_loop_stable(k2)
-
-    # rb=0 (on z-axis) requires special treatment because ellipe(x)/x and
-    # cel(x)/x for x->0 both appear in the expressions, both are finite at x=0
-    # but evaluation gives 0/0=nan
-    # To avoid treating x=0 as a special case (with masks), one would have to find
-    # designated algorithms for these expressions (which require additional evaluation
-    # and computation time because xi must be evaluated in any case).
-
-    mask1 = rb == 0
-    z_over_r = np.zeros(nX)
-    k2_over_rb = np.ones(nX) * 4 / (z2 + 1)
-    z_over_r[~mask1] = zb[~mask1] / rb[~mask1]  # will be zero when r=0
-    k2_over_rb[~mask1] = k2[~mask1] / rb[~mask1]  # will be zero when r=0
-
-    # field components
-    pf = 1 / np.sqrt(brack) / (1 - k2) / rad
-    Br = pf * z_over_r * xi
-    Bz = pf * (k2_over_rb * ellipe(k2) - xi)
-
-    # current and [mT] unit
+    r0 = np.abs(diameter / 2)
+    n = len(r0)
+
+    # allocate result
     Br_tot, Bz_tot = np.zeros((2, n))
-    Br_tot[mask_general] = Br
-    Bz_tot[mask_general] = Bz
+
+    # Special cases:
+    # case1: loop radius is 0 -> return (0,0,0)
+    mask1 = r0 == 0
+    # case2: at singularity -> return (0,0,0)
+    mask2 = np.logical_and(abs(r - r0) < 1e-15 * r0, z == 0)
+    # case3: r=0
+    mask3 = r == 0
+    if np.any(mask3):
+        mask4 = mask3 * ~mask1  # only relevant if not also case1
+        Bz_tot[mask4] = (
+            0.6283185307179587
+            * r0[mask4] ** 2
+            / (z[mask4] ** 2 + r0[mask4] ** 2) ** (3 / 2)
+        )
+
+    # general case
+    mask5 = ~np.logical_or(np.logical_or(mask1, mask2), mask3)
+    if np.any(mask5):
+        r0 = r0[mask5]
+        r = r[mask5]
+        z = z[mask5]
+        n5 = len(r0)
+
+        # express through ratios (make dimensionless, avoid large/small input values, stupid)
+        r = r / r0
+        z = z / r0
+
+        # field computation from paper
+        z2 = z**2
+        x0 = z2 + (r + 1) ** 2
+        k2 = 4 * r / x0
+        q2 = (z2 + (r - 1) ** 2) / x0
+
+        k = np.sqrt(k2)
+        q = np.sqrt(q2)
+        p = 1 + q
+        pf = k / np.sqrt(r) / q2 / 20 / r0
+
+        # cel* part
+        cc = k2 * k2
+        ss = 2 * cc * q / p
+        Br_tot[mask5] = pf * z / r * cel_iter(q, p, np.ones(n5), cc, ss, p, q)
+
+        # cel** part
+        cc = k2 * (k2 - (q2 + 1) / r)
+        ss = 2 * k2 * q * (k2 / p - p / r)
+        Bz_tot[mask5] = -pf * cel_iter(q, p, np.ones(n5), cc, ss, p, q)
 
     # transform field to cartesian CS
     Bx_tot, By_tot = cyl_field_to_cart(phi, Br_tot)
@@ -133,6 +127,6 @@ def current_loop_field(
 
     # B or H field
     if bh:
-        return B_cart / 10
+        return B_cart
 
-    return B_cart / 4 / np.pi
+    return B_cart / np.pi * 2.5
diff --git a/magpylib/_src/fields/field_BH_sphere.py b/magpylib/_src/fields/field_BH_sphere.py
index 1d1b5860a..6a30352ed 100644
--- a/magpylib/_src/fields/field_BH_sphere.py
+++ b/magpylib/_src/fields/field_BH_sphere.py
@@ -7,6 +7,7 @@
 from magpylib._src.input_checks import check_field_input
 
 
+# CORE
 def magnet_sphere_field(
     field: str,
     observers: np.ndarray,
diff --git a/magpylib/_src/fields/field_BH_tetrahedron.py b/magpylib/_src/fields/field_BH_tetrahedron.py
new file mode 100644
index 000000000..ec2a991c3
--- /dev/null
+++ b/magpylib/_src/fields/field_BH_tetrahedron.py
@@ -0,0 +1,145 @@
+"""
+Implementation for the magnetic field of homogeneously
+magnetized tetrahedra. Computation details in function docstrings.
+"""
+import numpy as np
+
+from magpylib._src.fields.field_BH_triangle import triangle_field
+from magpylib._src.input_checks import check_field_input
+
+
+def check_chirality(points: np.ndarray) -> np.ndarray:
+    """
+    Checks if quartupel of points (p0,p1,p2,p3) that forms tetrahedron is arranged in a way
+    that the vectors p0p1, p0p2, p0p3 form a right-handed system
+
+    Parameters
+    -----------
+    points: 3d-array of shape (m x 4 x 3)
+            m...number of tetrahedrons
+
+    Returns
+    ----------
+    new list of points, where p2 and p3 are possibly exchanged so that all
+    tetrahedron is given in a right-handed system.
+    """
+
+    vecs = np.zeros((len(points), 3, 3))
+    vecs[:, :, 0] = points[:, 1, :] - points[:, 0, :]
+    vecs[:, :, 1] = points[:, 2, :] - points[:, 0, :]
+    vecs[:, :, 2] = points[:, 3, :] - points[:, 0, :]
+
+    dets = np.linalg.det(vecs)
+    dets_neg = dets < 0
+
+    if np.any(dets_neg):
+        points[dets_neg, 2:, :] = points[dets_neg, 3:1:-1, :]
+
+    return points
+
+
+def point_inside(points: np.ndarray, vertices: np.ndarray) -> np.ndarray:
+    """
+    Takes points, as well as the vertices of a tetrahedra.
+    Returns boolean array indicating whether the points are inside the tetrahedra.
+    """
+    mat = vertices[:, 1:].swapaxes(0, 1) - vertices[:, 0]
+    mat = np.transpose(mat.swapaxes(0, 1), (0, 2, 1))
+
+    tetra = np.linalg.inv(mat)
+    newp = np.matmul(tetra, np.reshape(points - vertices[:, 0, :], (*points.shape, 1)))
+    inside = (
+        np.all(newp >= 0, axis=1)
+        & np.all(newp <= 1, axis=1)
+        & (np.sum(newp, axis=1) <= 1)
+    ).flatten()
+
+    return inside
+
+
+def magnet_tetrahedron_field(
+    field: str,
+    observers: np.ndarray,
+    magnetization: np.ndarray,
+    vertices: np.ndarray,
+) -> np.ndarray:
+    """
+    Magnetic field generated by a homogeneously magnetized tetrahedron.
+
+    Parameters
+    ----------
+    field: str, default=`'B'`
+        If `field='B'` return B-field in units of [mT], if `field='H'` return H-field
+        in units of [kA/m].
+
+    observers: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    magnetization: ndarray, shape (n,3)
+        Homogeneous magnetization vector in units of [mT].
+
+    vertices: ndarray, shape (n,4,3)
+        Vertices of the individual tetrahedrons [(pos1a, pos1b, pos1c, pos1d),
+        (pos2a, pos2b, pos2c, pos2d), ...].
+
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of magnet in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
+
+    Examples
+    --------
+    Compute the B-field of two different tetrahedron-observer instances
+
+    >>> import numpy as np
+    >>> import magpylib as magpy
+    >>>
+    >>> obs = np.array([(1,2,3), (2,3,4)])
+    >>> mag = np.array([(222,333,444), (111,112,113)])
+    >>> vert = np.array([((-1,0,0), (1,-1,0), (1,1,0), (0,0,1))]*2)
+    >>> B = magpy.core.magnet_tetrahedron_field('B', obs, mag, vert)
+    >>> print(B)
+    [[0.19075398 0.8240532  1.18170862]
+     [0.03125701 0.08445416 0.1178967 ]]
+
+    Notes
+    -----
+    The tetrahedron is built up via 4 faces applying the Triangle class, making sure that
+    all normal vectors point outwards, and providing inside-outside evaluation to distinguish
+    between B- and H-field.
+    """
+
+    bh = check_field_input(field, "magnet_tetrahedron_field()")
+
+    n = len(observers)
+
+    vertices = check_chirality(vertices)
+    tri_vertices = np.concatenate(
+        (
+            vertices[:, (0, 2, 1), :],
+            vertices[:, (0, 1, 3), :],
+            vertices[:, (1, 2, 3), :],
+            vertices[:, (0, 3, 2), :],
+        ),
+        axis=0,
+    )
+    tri_fields = triangle_field(
+        field,
+        np.tile(observers, (4, 1)),
+        np.tile(magnetization, (4, 1)),
+        tri_vertices,
+    )
+    tetra_field = (  # slightly faster than reshape + sum
+        tri_fields[:n]
+        + tri_fields[n : 2 * n]
+        + tri_fields[2 * n : 3 * n]
+        + tri_fields[3 * n :]
+    )
+
+    if not bh:
+        return tetra_field
+
+    # if B, and inside magnet add magnetization vector
+    mask_inside = point_inside(observers, vertices)
+    tetra_field[mask_inside] += magnetization[mask_inside]
+    return tetra_field
diff --git a/magpylib/_src/fields/field_BH_triangle.py b/magpylib/_src/fields/field_BH_triangle.py
new file mode 100644
index 000000000..213415144
--- /dev/null
+++ b/magpylib/_src/fields/field_BH_triangle.py
@@ -0,0 +1,186 @@
+"""
+Implementations of analytical expressions for the magnetic field of a triangular surface.
+Computation details in function docstrings.
+"""
+# pylance: disable=Code is unreachable
+import numpy as np
+
+from magpylib._src.input_checks import check_field_input
+
+
+def vcross3(a: np.ndarray, b: np.ndarray) -> np.ndarray:
+    """
+    vectorized cross product for 3d vectors. Is ~4x faster than np.cross when
+    arrays are smallish. Only slightly faster for large arrays.
+    input shape a,b: (n,3)
+    returns: (n, 3)
+    """
+    # receives nan values at corners
+    with np.errstate(invalid="ignore"):
+        result = np.array(
+            [
+                a[:, 1] * b[:, 2] - a[:, 2] * b[:, 1],
+                a[:, 2] * b[:, 0] - a[:, 0] * b[:, 2],
+                a[:, 0] * b[:, 1] - a[:, 1] * b[:, 0],
+            ]
+        )
+    return result.T
+
+
+def norm_vector(v) -> np.ndarray:
+    """
+    Calculates normalized orthogonal vector on a plane defined by three vertices.
+    """
+    a = v[:, 1] - v[:, 0]
+    b = v[:, 2] - v[:, 0]
+    n = vcross3(a, b)
+    n_norm = np.linalg.norm(n, axis=-1)
+    return n / np.expand_dims(n_norm, axis=-1)
+
+
+def solid_angle(R: np.ndarray, r: np.ndarray) -> np.ndarray:
+    """
+    Vectorized computation of the solid angle of triangles.
+
+    Triangle point indices are 1,2,3, different triangles are denoted by a,b,c,...
+    The first triangle is defined as R1a, R2a, R3a.
+
+    Input:
+    R = [(R1a, R1b, R1c, ...), (R2a, R2b, R2c, ...), (R3a, R3b, R3c, ...)]
+    r = [(|R1a|, |R1b|, |R1c|, ...), (|R2a|, |R2b|, |R2c|, ...), (|R3a|, |R3b|, |R3c|, ...)]
+
+    Returns:
+    [sangle_a, sangle_b, sangle_c, ...]
+    """
+
+    # Calculates (oriented) volume of the parallelepiped in vectorized form.
+    N = np.einsum("ij, ij->i", R[2], vcross3(R[1], R[0]))
+
+    D = (
+        r[0] * r[1] * r[2]
+        + np.einsum("ij, ij->i", R[2], R[1]) * r[0]
+        + np.einsum("ij, ij->i", R[2], R[0]) * r[1]
+        + np.einsum("ij, ij->i", R[1], R[0]) * r[2]
+    )
+    result = 2.0 * np.arctan2(N, D)
+
+    # modulus 2pi to avoid jumps on edges in line
+    # "B = sigma * ((n.T * solid_angle(R, r)) - vcross3(n, PQR).T)"
+    # <-- bad fix :(
+
+    return np.where(abs(result) > 6.2831853, 0, result)
+
+
+def triangle_field(
+    field: str,
+    observers: np.ndarray,
+    magnetization: np.ndarray,
+    vertices: np.ndarray,
+) -> np.ndarray:
+    """Magnetic field generated by homogeneous (magnetic) charges on triangular surfaces.
+
+    Can be used to compute the field of a homogeneously magnetized bodies
+    with triangular surface mesh. In this case each Triangle must be defined so that the
+    normal vector points outwards.
+
+    Parameters
+    ----------
+    field: str, default=`'B'`
+        If `field='B'` return B-field in units of [mT], if `field='H'` return H-field
+        in units of [kA/m].
+
+    observers: ndarray, shape (n,3)
+        Observer positions (x,y,z) in Cartesian coordinates in units of [mm].
+
+    magnetization: ndarray, shape (n,3)
+        Homogeneous magnetization vector in units of [mT]. The triangle surface charge is the
+        projection of the magnetization vector onto the surface normal vector. The order
+        of the vertices defines the sign of the normal vector (right-hand-rule).
+
+    vertices: ndarray, shape (n,3,3)
+        Vertices of triangular faces of the format
+        [(pos1a, pos1b, pos1c), (pos2a, pos2b, pos2c), ...].
+        The vertex order defines the normal vector orientation of the face (right-hand-rule).
+
+    Returns
+    -------
+    B-field or H-field: ndarray, shape (n,3)
+        B/H-field of magnet in Cartesian coordinates (Bx, By, Bz) in units of [mT]/[kA/m].
+
+    Examples
+    --------
+    Compute the B-field of two different triangle-observer instances.
+
+    >>> import numpy as np
+    >>> import magpylib as magpy
+
+    >>> mag = np.array([(22,33,44), (33,44,55)])
+    >>> vert = np.array([((-1,0,0), (1,-1,0), (1,1,0))]*2)
+    >>> obs = np.array([(1,2,3), (2,3,4)])
+    >>> B = magpy.core.triangle_field('B', obs, mag, vert)
+    >>> print(B)
+    [[0.08836711 0.27062047 0.42198603]
+     [0.09721743 0.17609525 0.23937461]]
+
+    Notes
+    -----
+    Field computations implemented from Guptasarma, Geophysics, 1999, 64:1, 70-74.
+    Corners give (nan, nan, nan). Edges and in-plane perp components are set to 0.
+    Loss of precision when approaching a triangle as (x-edge)**2 :(
+    Loss of precision with distance from the triangle as distance**3 :(
+    """
+    # pylint: disable=too-many-statements
+    bh = check_field_input(field, "triangle_field()")
+
+    n = norm_vector(vertices)
+    sigma = np.einsum("ij, ij->i", n, magnetization)  # vectorized inner product
+
+    # vertex <-> observer
+    R = np.swapaxes(vertices, 0, 1) - observers
+    r2 = np.sum(R * R, axis=-1)
+    r = np.sqrt(r2)
+
+    # vertex <-> vertex
+    L = vertices[:, (1, 2, 0)] - vertices[:, (0, 1, 2)]
+    L = np.swapaxes(L, 0, 1)
+    l2 = np.sum(L * L, axis=-1)
+    l = np.sqrt(l2)
+
+    # vert-vert -- vert-obs
+    b = np.einsum("ijk, ijk->ij", R, L)
+    bl = b / l
+    ind = np.fabs(r + bl)  # closeness measure to corner and edge
+
+    # The computation of ind is the origin of a major numerical instability
+    #    when approaching the triangle because r ~ -bl. This number
+    #    becomes small at the same rate as it looses precision.
+    #    This is a major problem, because at distances 1e-8 and 1e8 all precision
+    #    is already lost !!!
+    # The second problem is at corner and edge extensions where ind also computes
+    #    as 0. Here one approaches a special case where another evaluation should
+    #    be used. This problem is solved in the following lines.
+    # np.seterr must be used because of a numpy bug. It does not interpret where
+    #   correctly. The following code will raise a numpy warning - but obviously shouldn't
+    #
+    # x = np.array([(0,1,2), (0,0,1)])
+    # np.where(
+    #     x>0,
+    #     1/x,
+    #     0
+    # )
+
+    with np.errstate(divide="ignore", invalid="ignore"):
+        I = np.where(
+            ind > 1.0e-12,
+            1.0 / l * np.log((np.sqrt(l2 + 2 * b + r2) + l + bl) / ind),
+            -(1.0 / l) * np.log(np.fabs(l - r) / r),
+        )
+    PQR = np.einsum("ij, ijk -> jk", I, L)
+    B = sigma * (n.T * solid_angle(R, r) - vcross3(n, PQR).T)
+
+    # return B or compute and return H -------------
+    if bh:
+        return B.T / np.pi / 4.0
+
+    H = B.T / 1.6 / np.pi**2  # mT -> kA/m
+    return H
diff --git a/magpylib/_src/fields/special_cel.py b/magpylib/_src/fields/special_cel.py
index c144f046a..f9690fd4a 100644
--- a/magpylib/_src/fields/special_cel.py
+++ b/magpylib/_src/fields/special_cel.py
@@ -1,3 +1,5 @@
+import math as m
+
 import numpy as np
 
 
@@ -130,73 +132,53 @@ def ellipticPi(x, y):
     return celv(kcv, pv, cv, sv)
 
 
-def cel_loop_stable(k2):
+def cel_iter(qc, p, g, cc, ss, em, kk):
     """
-    numerically stabilized version of the function
-    xi_loop = (2-k**2)E(k**2)-2(1-k**2)K(k**2) = cel(np.sqrt(1-k2), 1, k2, k2*(k2-1))
-    needed for the circular current loop field. This modicfication of the cel()
-    algorithm is numerically stable when k2 becomes small.
+    Iterative part of Bulirsch cel algorithm
     """
-    # if k2 == 0: return 0    # on axis
+    # case1: scalar input
+    #   This cannot happen in core functions
 
-    pp = 1 - k2  # allocate pp and use temporarily for 1-k2
-    k = np.sqrt(pp)
-    kk = 1 + k  # allocate pp and use temporarily for 1+k
+    # case2: small input vector - loop is faster than vectorized computation
+    n_input = len(qc)
+    if n_input < 15:
+        result = np.zeros(n_input)
+        for i in range(n_input):
+            result[i] = cel_iter0(qc[i], p[i], g[i], cc[i], ss[i], em[i], kk[i])
 
-    g = 1 if isinstance(k2, (float, int)) else np.ones(len(k2))
-    cc = k2**2
-    ss = 2 * cc * pp / (kk - k2)
-    pp = kk
-    em = kk
-    kk = k
+    # case3: vectorized evaluation
+    return cel_iterv(qc, p, g, cc, ss, em, kk)
 
-    errtol = 0.000001
-    while np.any(abs(g - k) > g * errtol):
-        k = 2 * np.sqrt(kk)
-        kk = k * em
+
+def cel_iter0(qc, p, g, cc, ss, em, kk):
+    """
+    Iterative part of Bulirsch cel algorithm
+    """
+    while m.fabs(g - qc) >= qc * 1e-8:
+        qc = 2 * m.sqrt(kk)
+        kk = qc * em
         f = cc
-        cc = cc + ss / pp
-        g = kk / pp
+        cc = cc + ss / p
+        g = kk / p
         ss = 2 * (ss + f * g)
-        pp = g + pp
+        p = p + g
         g = em
-        em = k + em
-    return (np.pi / 2) * (ss + cc * em) / (em * (em + pp))
+        em = em + qc
+    return 1.5707963267948966 * (ss + cc * em) / (em * (em + p))
 
 
-# def cel_loop_stable_old(k2):
-#     """
-#     numerically stabilized version of the function
-#         xi_loop = (2-k**2)E(k**2)-2(1-k**2)K(k**2)
-#         needed for the current.loop fields. See paper
-#         Leitner2021
-#     """
-#     n = len(k2)
-#     result = np.empty(n)
-
-#     mask1 = (k2 > 0.04)
-
-#     if np.any(mask1):
-#         k2m = k2[mask1]
-#         result[mask1] = cel(np.sqrt(1-k2m), np.ones(np.sum(mask1)), k2m, k2m*(k2m-1))
-#     if np.any(~mask1):
-#         result[~mask1] = cel_loop_taylor0(k2[~mask1])
-#     return result
-
-
-# def cel_loop_taylor0(k2):
-#     """
-#     taylor expansion of the function xi_loop about k2=0
-#         See paper Leitner2021
-#     """
-#     C2 = 0.5890486225480862
-#     C3 = 0.1472621556370216
-#     C4 = 0.06902913545485386
-#     C5 = 0.04026699568199808
-#     C6 = 0.02642521591631124
-#     C7 = 0.01868640268367724
-#     C8 = 0.01391747699878044
-#     C9 = 0.01076947624905629
-#     C10 = 0.008581926385966734
-#     return (C2*k2**2 + C3*k2**3 + C4*k2**4 + C5*k2**5
-#         + C6*k2**6 + C7*k2**7 + C8*k2**8 + C9*k2**9 + C10*k2**10)
+def cel_iterv(qc, p, g, cc, ss, em, kk):
+    """
+    Iterative part of Bulirsch cel algorithm
+    """
+    while np.any(np.fabs(g - qc) >= qc * 1e-8):
+        qc = 2 * np.sqrt(kk)
+        kk = qc * em
+        f = cc
+        cc = cc + ss / p
+        g = kk / p
+        ss = 2 * (ss + f * g)
+        p = p + g
+        g = em
+        em = em + qc
+    return 1.5707963267948966 * (ss + cc * em) / (em * (em + p))
diff --git a/magpylib/_src/input_checks.py b/magpylib/_src/input_checks.py
index e139f6ebc..4704b6f90 100644
--- a/magpylib/_src/input_checks.py
+++ b/magpylib/_src/input_checks.py
@@ -47,7 +47,7 @@ def make_float_array(inp, msg: str):
     return inp_array
 
 
-def check_array_shape(inp: np.ndarray, dims: tuple, shape_m1: int, msg: str):
+def check_array_shape(inp: np.ndarray, dims: tuple, shape_m1: int, length=None, msg=""):
     """check if inp shape is allowed
     inp: test object
     dims: list, list of allowed dims
@@ -55,9 +55,12 @@ def check_array_shape(inp: np.ndarray, dims: tuple, shape_m1: int, msg: str):
     msg: str, error msg
     """
     if inp.ndim in dims:
-        if inp.shape[-1] == shape_m1:
-            return None
-        if shape_m1 == "any":
+        if length is None:
+            if inp.shape[-1] == shape_m1:
+                return None
+            if shape_m1 == "any":
+                return None
+        elif len(inp) == length:
             return None
     raise MagpylibBadUserInput(msg)
 
@@ -309,6 +312,7 @@ def check_format_input_vector(
     shape_m1,
     sig_name,
     sig_type,
+    length=None,
     reshape=False,
     allow_None=False,
     forbid_negative0=False,
@@ -339,13 +343,15 @@ def check_format_input_vector(
         inp,
         dims=dims,
         shape_m1=shape_m1,
+        length=length,
         msg=(
             f"Input parameter `{sig_name}` must be {sig_type}.\n"
             f"Instead received array_like with shape {inp.shape}."
         ),
     )
-    if reshape:
-        return np.reshape(inp, (-1, 3))
+    if isinstance(reshape, tuple):
+        return np.reshape(inp, reshape)
+
     if forbid_negative0:
         if np.any(inp <= 0):
             raise MagpylibBadUserInput(
@@ -555,6 +561,7 @@ def check_dimensions(sources):
     from magpylib._src.obj_classes.class_magnet_Sphere import Sphere
     from magpylib._src.obj_classes.class_current_Line import Line
     from magpylib._src.obj_classes.class_current_Loop import Loop
+    from magpylib._src.obj_classes.class_magnet_Tetrahedron import Tetrahedron
 
     for src in sources:
         if isinstance(src, (Cuboid, Cylinder, CylinderSegment)):
@@ -567,7 +574,7 @@ def check_dimensions(sources):
                 raise MagpylibMissingInput(
                     f"Parameter `diameter` of {src} must be set."
                 )
-        elif isinstance(src, Line):
+        elif isinstance(src, (Line, Tetrahedron)):
             if src.vertices is None:
                 raise MagpylibMissingInput(
                     f"Parameter `vertices` of {src} must be set."
diff --git a/magpylib/_src/obj_classes/class_BaseDisplayRepr.py b/magpylib/_src/obj_classes/class_BaseDisplayRepr.py
index 71b65070e..3c647512f 100644
--- a/magpylib/_src/obj_classes/class_BaseDisplayRepr.py
+++ b/magpylib/_src/obj_classes/class_BaseDisplayRepr.py
@@ -1,6 +1,4 @@
-"""BaseGeo class code
-READY FOR V4
-"""
+"""BaseGeo class code"""
 # pylint: disable=cyclic-import
 import numpy as np
 
diff --git a/magpylib/_src/obj_classes/class_BaseExcitations.py b/magpylib/_src/obj_classes/class_BaseExcitations.py
index 374adfa12..4b0f674c0 100644
--- a/magpylib/_src/obj_classes/class_BaseExcitations.py
+++ b/magpylib/_src/obj_classes/class_BaseExcitations.py
@@ -1,6 +1,4 @@
-"""BaseHomMag class code
-DOCSTRINGS V4 READY
-"""
+"""BaseHomMag class code"""
 # pylint: disable=cyclic-import
 from magpylib._src.fields.field_wrap_BH import getBH_level2
 from magpylib._src.input_checks import check_format_input_scalar
diff --git a/magpylib/_src/obj_classes/class_BaseGeo.py b/magpylib/_src/obj_classes/class_BaseGeo.py
index 7a7b2e619..4ccced23e 100644
--- a/magpylib/_src/obj_classes/class_BaseGeo.py
+++ b/magpylib/_src/obj_classes/class_BaseGeo.py
@@ -1,6 +1,4 @@
-"""BaseGeo class code
-DOCSTRING v4 READY
-"""
+"""BaseGeo class code"""
 # pylint: disable=cyclic-import
 # pylint: disable=too-many-instance-attributes
 # pylint: disable=protected-access
@@ -113,7 +111,7 @@ def _init_position_orientation(self, position, orientation):
             shape_m1=3,
             sig_name="position",
             sig_type="array_like (list, tuple, ndarray) with shape (3,) or (n,3)",
-            reshape=True,
+            reshape=(-1, 3),
         )
         oriQ = check_format_input_orientation(orientation, init_format=True)
 
@@ -183,7 +181,7 @@ def position(self, inp):
             shape_m1=3,
             sig_name="position",
             sig_type="array_like (list, tuple, ndarray) with shape (3,) or (n,3)",
-            reshape=True,
+            reshape=(-1, 3),
         )
 
         # pad/slice and set orientation path to same length
diff --git a/magpylib/_src/obj_classes/class_BaseTransform.py b/magpylib/_src/obj_classes/class_BaseTransform.py
index 62626ee87..d1eab5097 100644
--- a/magpylib/_src/obj_classes/class_BaseTransform.py
+++ b/magpylib/_src/obj_classes/class_BaseTransform.py
@@ -1,6 +1,4 @@
-"""BaseTransform class code
-DOCSTRING v4 READY
-"""
+"""BaseTransform class code"""
 # pylint: disable=protected-access
 import numbers
 
diff --git a/magpylib/_src/obj_classes/class_Collection.py b/magpylib/_src/obj_classes/class_Collection.py
index df06ad538..acc21845a 100644
--- a/magpylib/_src/obj_classes/class_Collection.py
+++ b/magpylib/_src/obj_classes/class_Collection.py
@@ -739,8 +739,8 @@ class Collection(BaseGeo, BaseCollection):
 
     The field can be computed at position (0,0,0) as if the collection was a single source:
 
-    B = col.getB((0,0,0))
-    print(B)
+    >>> B = col.getB((0,0,0))
+    >>> print(B)
     [ 0.00126232 -0.00093169 -0.00034448]
 
     We add a sensor at position (0,0,0) to the collection:
@@ -754,8 +754,8 @@ class Collection(BaseGeo, BaseCollection):
     and can compute the field of the sources in the collection seen by the sensor with
     a single command:
 
-    B = col.getB()
-    print(B)
+    >>> B = col.getB()
+    >>> print(B)
     [ 0.00126232 -0.00093169 -0.00034448]
     """
 
diff --git a/magpylib/_src/obj_classes/class_Sensor.py b/magpylib/_src/obj_classes/class_Sensor.py
index 7ead0500e..7831a741b 100644
--- a/magpylib/_src/obj_classes/class_Sensor.py
+++ b/magpylib/_src/obj_classes/class_Sensor.py
@@ -1,6 +1,4 @@
-"""Sensor class code
-DOCSTRINGS V4 READY
-"""
+"""Sensor class code"""
 from magpylib._src.display.traces_generic import make_Sensor
 from magpylib._src.fields.field_wrap_BH import getBH_level2
 from magpylib._src.input_checks import check_format_input_vector
diff --git a/magpylib/_src/obj_classes/class_current_Line.py b/magpylib/_src/obj_classes/class_current_Line.py
index d7d66fb2d..bc5859081 100644
--- a/magpylib/_src/obj_classes/class_current_Line.py
+++ b/magpylib/_src/obj_classes/class_current_Line.py
@@ -1,6 +1,4 @@
-"""Line current class code
-DOCSTRINGS V4 READY
-"""
+"""Line current class code"""
 from magpylib._src.display.traces_generic import make_Line
 from magpylib._src.fields.field_BH_line import current_vertices_field
 from magpylib._src.input_checks import check_format_input_vertices
@@ -63,7 +61,7 @@ class Line(BaseCurrent):
     We rotate the source object, and compute the B-field, this time at a set of observer positions:
 
     >>> src.rotate_from_angax(45, 'x')
-    Line...
+    Line(id=...)
     >>> B = src.getB([(1,1,1), (2,2,2), (3,3,3)])
     >>> print(B)
     [[-6.68990257e-18  3.50341393e-02 -3.50341393e-02]
@@ -74,7 +72,7 @@ class Line(BaseCurrent):
     observer at position (1,1,1). This time we use a `Sensor` object as observer.
 
     >>> src.move([(-1,-1,-1), (-2,-2,-2)])
-    Line...
+    Line(id=...)
     >>> sens = magpy.Sensor(position=(1,1,1))
     >>> B = src.getB(sens)
     >>> print(B)
diff --git a/magpylib/_src/obj_classes/class_current_Loop.py b/magpylib/_src/obj_classes/class_current_Loop.py
index d63b7ee5c..f883ec9da 100644
--- a/magpylib/_src/obj_classes/class_current_Loop.py
+++ b/magpylib/_src/obj_classes/class_current_Loop.py
@@ -1,6 +1,4 @@
-"""Loop current class code
-DOCSTRINGS V4 READY
-"""
+"""Loop current class code"""
 from magpylib._src.display.traces_generic import make_Loop
 from magpylib._src.fields.field_BH_loop import current_loop_field
 from magpylib._src.input_checks import check_format_input_scalar
diff --git a/magpylib/_src/obj_classes/class_magnet_Cuboid.py b/magpylib/_src/obj_classes/class_magnet_Cuboid.py
index 01d66ba10..7b86b91a8 100644
--- a/magpylib/_src/obj_classes/class_magnet_Cuboid.py
+++ b/magpylib/_src/obj_classes/class_magnet_Cuboid.py
@@ -1,6 +1,4 @@
-"""Magnet Cuboid class code
-DOCSTRINGS V4 READY
-"""
+"""Magnet Cuboid class code"""
 from magpylib._src.display.traces_generic import make_Cuboid
 from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
 from magpylib._src.input_checks import check_format_input_vector
diff --git a/magpylib/_src/obj_classes/class_magnet_Cylinder.py b/magpylib/_src/obj_classes/class_magnet_Cylinder.py
index a42208cdb..7d5907a6a 100644
--- a/magpylib/_src/obj_classes/class_magnet_Cylinder.py
+++ b/magpylib/_src/obj_classes/class_magnet_Cylinder.py
@@ -1,6 +1,4 @@
-"""Magnet Cylinder class code
-DOCSTRINGS V4 READY
-"""
+"""Magnet Cylinder class code"""
 from magpylib._src.display.traces_generic import make_Cylinder
 from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_field
 from magpylib._src.input_checks import check_format_input_vector
diff --git a/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py b/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py
index 255aa31c2..62b60ab94 100644
--- a/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py
+++ b/magpylib/_src/obj_classes/class_magnet_CylinderSegment.py
@@ -1,6 +1,4 @@
-"""Magnet Cylinder class code
-DOCSTRINGS V4 READY
-"""
+"""Magnet Cylinder class code"""
 import numpy as np
 
 from magpylib._src.display.traces_generic import make_CylinderSegment
@@ -37,6 +35,10 @@ class CylinderSegment(BaseMagnet):
         Object position(s) in the global coordinates in units of [mm]. For m>1, the
         `position` and `orientation` attributes together represent an object path.
 
+    barycenter: array_like, shape (3,)
+        Read only property that returns the geometric barycenter (=center of mass)
+        of the object.
+
     orientation: scipy `Rotation` object with length 1 or m, default=`None`
         Object orientation(s) in the global coordinates. `None` corresponds to
         a unit-rotation. For m>1, the `position` and `orientation` attributes
diff --git a/magpylib/_src/obj_classes/class_magnet_Sphere.py b/magpylib/_src/obj_classes/class_magnet_Sphere.py
index 154f837a2..43ca2845c 100644
--- a/magpylib/_src/obj_classes/class_magnet_Sphere.py
+++ b/magpylib/_src/obj_classes/class_magnet_Sphere.py
@@ -1,6 +1,4 @@
-"""Magnet Sphere class code
-DOCSTRINGS V4 READY
-"""
+"""Magnet Sphere class code"""
 from magpylib._src.display.traces_generic import make_Sphere
 from magpylib._src.fields.field_BH_sphere import magnet_sphere_field
 from magpylib._src.input_checks import check_format_input_scalar
diff --git a/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py b/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py
new file mode 100644
index 000000000..e76cf6ae6
--- /dev/null
+++ b/magpylib/_src/obj_classes/class_magnet_Tetrahedron.py
@@ -0,0 +1,147 @@
+"""Magnet Tetrahedron class code"""
+import numpy as np
+
+from magpylib._src.display.traces_generic import make_Tetrahedron
+from magpylib._src.fields.field_BH_tetrahedron import magnet_tetrahedron_field
+from magpylib._src.input_checks import check_format_input_vector
+from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
+
+
+class Tetrahedron(BaseMagnet):
+    """Tetrahedron magnet with homogeneous magnetization.
+
+    Can be used as `sources` input for magnetic field computation.
+
+    When `position=(0,0,0)` and `orientation=None` the Tetrahedron vertices coordinates
+    are the same as in the global coordinate system. The geometric center of the Tetrahedron
+    is determined by its vertices and. It is not necessarily located in the origin an can
+    be computed with the barycenter property.
+
+    Parameters
+    ----------
+    magnetization: array_like, shape (3,), default=`None`
+        Magnetization vector (mu0*M, remanence field) in units of [mT] given in
+        the local object coordinates (rotates with object).
+
+    vertices: ndarray, shape (4,3)
+        Vertices [(x1,y1,z1), (x2,y2,z2), (x3,y3,z3), (x4,y4,z4)], in the relative
+        coordinate system of the tetrahedron.
+
+    position: array_like, shape (3,) or (m,3)
+        Object position(s) in the global coordinates in units of [mm]. For m>1, the
+        `position` and `orientation` attributes together represent an object path.
+        When setting vertices, the initial position is set to the barycenter.
+
+    barycenter: array_like, shape (3,)
+        Read only property that returns the geometric barycenter (=center of mass)
+        of the object.
+
+    orientation: scipy `Rotation` object with length 1 or m, default=`None`
+        Object orientation(s) in the global coordinates. `None` corresponds to
+        a unit-rotation. For m>1, the `position` and `orientation` attributes
+        together represent an object path.
+
+    parent: `Collection` object or `None`
+        The object is a child of it's parent collection.
+
+    style: dict
+        Object style inputs must be in dictionary form, e.g. `{'color':'red'}` or
+        using style underscore magic, e.g. `style_color='red'`.
+
+    Returns
+    -------
+    magnet source: `Tetrahedron` object
+
+    Examples
+    --------
+    `Tetrahedron` magnets are magnetic field sources. Below we compute the H-field [kA/m] of a
+    tetrahedron magnet with magnetization (100,200,300) in units of [mT] dimensions defined
+    through the vertices (0,0,0), (1,0,0), (0,1,0) and (0,0,1) in units of [mm] at the
+    observer position (1,1,1) given in units of [mm]:
+
+    >>> import magpylib as magpy
+    >>> verts = [(0,0,0), (1,0,0), (0,1,0), (0,0,1)]
+    >>> src = magpy.magnet.Tetrahedron(magnetization=(100,200,300), vertices=verts)
+    >>> H = src.getH((1,1,1))
+    >>> print(H)
+    [2.07089783 1.65671826 1.2425387 ]
+
+    We rotate the source object, and compute the B-field, this time at a set of observer positions:
+
+    src.rotate_from_angax(45, 'x')
+    Tetrahedron(id=...)
+    B = src.getB([(1,1,1), (2,2,2), (3,3,3)])
+    print(B)
+    [[ 8.68006559e-01  2.00895792e+00 -5.03469140e-01]
+    [ 1.01357229e-01  1.93731796e-01 -1.59677364e-02]
+    [ 2.90426931e-02  5.22556994e-02 -1.70596096e-03]]
+
+    The same result is obtained when the rotated source moves along a path away from an
+    observer at position (1,1,1). Here we use a `Sensor` object as observer.
+
+    sens = magpy.Sensor(position=(1,1,1))
+    src.move([(-1,-1,-1), (-2,-2,-2)])
+    Sensor(id=...)
+    B = src.getB(sens)
+    print(B)
+    [[ 8.68006559e-01  2.00895792e+00 -5.03469140e-01]
+    [ 1.01357229e-01  1.93731796e-01 -1.59677364e-02]
+    [ 2.90426931e-02  5.22556994e-02 -1.70596096e-03]]
+    """
+
+    _field_func = staticmethod(magnet_tetrahedron_field)
+    _field_func_kwargs_ndim = {"magnetization": 1, "vertices": 3}
+    _draw_func = make_Tetrahedron
+
+    def __init__(
+        self,
+        magnetization=None,
+        vertices=None,
+        position=(0, 0, 0),
+        orientation=None,
+        style=None,
+        **kwargs,
+    ):
+
+        # instance attributes
+        self.vertices = vertices
+        self._object_type = "Tetrahedron"
+
+        # init inheritance
+        super().__init__(position, orientation, magnetization, style, **kwargs)
+
+    # property getters and setters
+    @property
+    def vertices(self):
+        """Length of the Tetrahedron sides [a,b,c] in units of [mm]."""
+        return self._vertices
+
+    @vertices.setter
+    def vertices(self, dim):
+        """Set Tetrahedron vertices (a,b,c), shape (3,), [mm]."""
+        self._vertices = check_format_input_vector(
+            dim,
+            dims=(2,),
+            shape_m1=3,
+            length=4,
+            sig_name="Tetrahedron.vertices",
+            sig_type="array_like (list, tuple, ndarray) of shape (4,3)",
+            allow_None=True,
+        )
+
+    @property
+    def _barycenter(self):
+        """Object barycenter."""
+        return self._get_barycenter(self._position, self._orientation, self.vertices)
+
+    @property
+    def barycenter(self):
+        """Object barycenter."""
+        return np.squeeze(self._barycenter)
+
+    @staticmethod
+    def _get_barycenter(position, orientation, vertices):
+        """Returns the barycenter of a tetrahedron."""
+        centroid = np.mean(vertices, axis=0)
+        barycenter = orientation.apply(centroid) + position
+        return barycenter
diff --git a/magpylib/_src/obj_classes/class_misc_Dipole.py b/magpylib/_src/obj_classes/class_misc_Dipole.py
index 15e96497d..cbebee5b5 100644
--- a/magpylib/_src/obj_classes/class_misc_Dipole.py
+++ b/magpylib/_src/obj_classes/class_misc_Dipole.py
@@ -1,6 +1,4 @@
-"""Dipole class code
-DOCSTRINGS V4 READY
-"""
+"""Dipole class code"""
 from magpylib._src.display.traces_generic import make_Dipole
 from magpylib._src.fields.field_BH_dipole import dipole_field
 from magpylib._src.input_checks import check_format_input_vector
diff --git a/magpylib/_src/obj_classes/class_misc_Triangle.py b/magpylib/_src/obj_classes/class_misc_Triangle.py
new file mode 100644
index 000000000..69d91077a
--- /dev/null
+++ b/magpylib/_src/obj_classes/class_misc_Triangle.py
@@ -0,0 +1,145 @@
+"""Magnet Triangle class
+"""
+import numpy as np
+
+from magpylib._src.display.traces_generic import make_Triangle
+from magpylib._src.fields.field_BH_triangle import triangle_field
+from magpylib._src.input_checks import check_format_input_vector
+from magpylib._src.obj_classes.class_BaseExcitations import BaseMagnet
+
+
+class Triangle(BaseMagnet):
+    """Triangular surface with homogeneous magnetic surface charge.
+
+    Can be used as `sources` input for magnetic field computation.
+
+    When `position=(0,0,0)` and `orientation=None` the local object coordinates of the
+    Triangle vertices coincide with the global coordinate system. The geometric
+    center of the Triangle is determined by its vertices.
+
+    Parameters
+    ----------
+    magnetization: array_like, shape (3,), default=`None`
+        Magnetization vector (mu0*M, remanence field) in units of [mT] given in
+        the local object coordinates (rotates with object). The homogeneous surface
+        charge of the Triangle is given by the projection of the magnetization on the
+        Triangle normal vector (right-hand-rule).
+
+    vertices: ndarray, shape (3,3)
+        Triple of vertices in the local object coordinates.
+
+    position: array_like, shape (3,) or (m,3)
+        Object position(s) in the global coordinates in units of [mm]. For m>1, the
+        `position` and `orientation` attributes together represent an object path.
+
+    barycenter: array_like, shape (3,)
+        Read only property that returns the geometric barycenter (=center of mass)
+        of the object.
+
+    orientation: scipy `Rotation` object with length 1 or m, default=`None`
+        Object orientation(s) in the global coordinates. `None` corresponds to
+        a unit-rotation. For m>1, the `position` and `orientation` attributes
+        together represent an object path.
+
+    parent: `Collection` object or `None`
+        The object is a child of it's parent collection.
+
+    style: dict
+        Object style inputs must be in dictionary form, e.g. `{'color':'red'}` or
+        using style underscore magic, e.g. `style_color='red'`.
+
+    Returns
+    -------
+    magnet source: `Triangle` object
+
+    Examples
+    --------
+    `Triangle` objects are magnetic field sources. Below we compute the H-field [kA/m] of a
+    Triangle object with magnetization (100,200,300) in units of [mT] dimensions defined
+    through the vertices (0,0,0), (1,0,0) and (0,1,0) in units of [mm] at the
+    observer position (1,1,1) given in units of [mm]:
+
+    >>> import magpylib as magpy
+    >>> verts = [(0,0,0), (1,0,0), (0,1,0)]
+    >>> src = magpy.misc.Triangle(magnetization=(100,200,300), vertices=verts)
+    >>> H = src.getH((1,1,1))
+    >>> print(H)
+    [2.24851814 2.24851814 3.49105683]
+
+    We rotate the source object, and compute the B-field, this time at a set of observer positions:
+
+    >>> src.rotate_from_angax(45, 'x')
+    Triangle(id=...)
+    >>> B = src.getB([(1,1,1), (2,2,2), (3,3,3)])
+    >>> print(B)
+    [[3.94659011 4.21772671 4.21772671]
+     [0.73745776 0.77325648 0.77325648]
+     [0.30049474 0.31043974 0.31043974]]
+
+    The same result is obtained when the rotated source moves along a path away from an
+    observer at position (1,1,1). Here we use a `Sensor` object as observer.
+
+    >>> sens = magpy.Sensor(position=(1,1,1))
+    >>> src.move([(-1,-1,-1), (-2,-2,-2)])
+    Triangle(id=...)
+    >>> B = src.getB(sens)
+    >>> print(B)
+    [[3.94659011 4.21772671 4.21772671]
+     [0.73745776 0.77325648 0.77325648]
+     [0.30049474 0.31043974 0.31043974]]
+
+    """
+
+    _field_func = staticmethod(triangle_field)
+    _field_func_kwargs_ndim = {"magnetization": 2, "vertices": 2}
+    _draw_func = make_Triangle
+
+    def __init__(
+        self,
+        magnetization=None,
+        vertices=None,
+        position=(0, 0, 0),
+        orientation=None,
+        style=None,
+        **kwargs,
+    ):
+
+        self.vertices = vertices
+
+        # init inheritance
+        super().__init__(position, orientation, magnetization, style, **kwargs)
+
+    # property getters and setters
+    @property
+    def vertices(self):
+        """Object faces"""
+        return self._vertices
+
+    @vertices.setter
+    def vertices(self, val):
+        """Set face vertices (a,b,c), shape (3,3), [mm]."""
+        self._vertices = check_format_input_vector(
+            val,
+            dims=(2,),
+            shape_m1=3,
+            sig_name="Triangle.vertices",
+            sig_type="array_like (list, tuple, ndarray) of shape (3,3)",
+            allow_None=True,
+        )
+
+    @property
+    def _barycenter(self):
+        """Object barycenter."""
+        return self._get_barycenter(self._position, self._orientation, self._vertices)
+
+    @property
+    def barycenter(self):
+        """Object barycenter."""
+        return np.squeeze(self._barycenter)
+
+    @staticmethod
+    def _get_barycenter(position, orientation, vertices):
+        """Returns the barycenter of the Triangle object."""
+        centroid = np.mean(vertices, axis=0)
+        barycenter = orientation.apply(centroid) + position
+        return barycenter
diff --git a/magpylib/core/__init__.py b/magpylib/core/__init__.py
index b19b728d9..bb4ce8359 100644
--- a/magpylib/core/__init__.py
+++ b/magpylib/core/__init__.py
@@ -10,6 +10,8 @@
     "magnet_cuboid_field",
     "magnet_cylinder_field",
     "magnet_cylinder_segment_field",
+    "triangle_field",
+    "magnet_tetrahedron_field",
 ]
 
 from magpylib._src.fields.field_BH_dipole import dipole_field
@@ -19,3 +21,5 @@
 from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
 from magpylib._src.fields.field_BH_cylinder import magnet_cylinder_field
 from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_field
+from magpylib._src.fields.field_BH_triangle import triangle_field
+from magpylib._src.fields.field_BH_tetrahedron import magnet_tetrahedron_field
diff --git a/magpylib/graphics/model3d/__init__.py b/magpylib/graphics/model3d/__init__.py
index f6ad0f308..7d315381a 100644
--- a/magpylib/graphics/model3d/__init__.py
+++ b/magpylib/graphics/model3d/__init__.py
@@ -11,6 +11,8 @@
     "make_Cuboid",
     "make_CylinderSegment",
     "make_Prism",
+    "make_Tetrahedron",
+    "make_TriangularMesh",
 ]
 
 from magpylib._src.display.traces_base import (
@@ -20,4 +22,6 @@
     make_Cuboid,
     make_CylinderSegment,
     make_Prism,
+    make_Tetrahedron,
+    make_TriangularMesh,
 )
diff --git a/magpylib/magnet/__init__.py b/magpylib/magnet/__init__.py
index bfa9dc7a7..407b6a015 100644
--- a/magpylib/magnet/__init__.py
+++ b/magpylib/magnet/__init__.py
@@ -2,10 +2,16 @@
 The `magpylib.magnet` subpackage contains all magnet classes.
 """
 
-__all__ = ["Cuboid", "Cylinder", "Sphere", "CylinderSegment"]
-
+__all__ = [
+    "Cuboid",
+    "Cylinder",
+    "Sphere",
+    "CylinderSegment",
+    "Tetrahedron",
+]
 
 from magpylib._src.obj_classes.class_magnet_Cuboid import Cuboid
 from magpylib._src.obj_classes.class_magnet_Cylinder import Cylinder
 from magpylib._src.obj_classes.class_magnet_Sphere import Sphere
 from magpylib._src.obj_classes.class_magnet_CylinderSegment import CylinderSegment
+from magpylib._src.obj_classes.class_magnet_Tetrahedron import Tetrahedron
diff --git a/magpylib/misc/__init__.py b/magpylib/misc/__init__.py
index 2ab07f329..5898a27fa 100644
--- a/magpylib/misc/__init__.py
+++ b/magpylib/misc/__init__.py
@@ -2,7 +2,8 @@
 The `magpylib.misc` sub-package contains miscellaneous source objects.
 """
 
-__all__ = ["Dipole", "CustomSource"]
+__all__ = ["Dipole", "CustomSource", "Triangle"]
 
 from magpylib._src.obj_classes.class_misc_Dipole import Dipole
 from magpylib._src.obj_classes.class_misc_CustomSource import CustomSource
+from magpylib._src.obj_classes.class_misc_Triangle import Triangle
diff --git a/setup.py b/setup.py
index a0bc24452..aacb66518 100644
--- a/setup.py
+++ b/setup.py
@@ -21,7 +21,7 @@
 import setuptools
 from setuptools.command.install import install
 
-_magPyVersion = "4.1.2"
+_magPyVersion = "4.2.0"
 _SphinxVersion = "4.4.0"
 _name = "magpylib"
 _description = "Free Python3 package to compute magnetic fields."


Test patch:
diff --git a/tests/test__missing_optional_modules.py b/tests/test__missing_optional_modules.py
index 4fc8bc005..7ee84dedf 100644
--- a/tests/test__missing_optional_modules.py
+++ b/tests/test__missing_optional_modules.py
@@ -1,5 +1,5 @@
 import sys
-from unittest import mock
+from unittest.mock import patch
 
 import pytest
 
@@ -9,7 +9,7 @@
 def test_show_with_missing_pyvista():
     """Should raise if pyvista is not installed"""
     src = magpy.magnet.Cuboid((0, 0, 1000), (1, 1, 1))
-    with mock.patch.dict(sys.modules, {"pyvista": None}):
+    with patch.dict(sys.modules, {"pyvista": None}):
         # with pytest.raises(ModuleNotFoundError):
         src.show(return_fig=True, backend="pyvista")
 
@@ -17,6 +17,6 @@ def test_show_with_missing_pyvista():
 def test_dataframe_output_missing_pandas():
     """test if pandas is installed when using dataframe output in `getBH`"""
     src = magpy.magnet.Cuboid((0, 0, 1000), (1, 1, 1))
-    with mock.patch.dict(sys.modules, {"pandas": None}):
+    with patch.dict(sys.modules, {"pandas": None}):
         with pytest.raises(ModuleNotFoundError):
             src.getB((0, 0, 0), output="dataframe")
diff --git a/tests/test_field_functions.py b/tests/test_core_field_functions.py
similarity index 60%
rename from tests/test_field_functions.py
rename to tests/test_core_field_functions.py
index 21af4dcdf..0e829f822 100644
--- a/tests/test_field_functions.py
+++ b/tests/test_core_field_functions.py
@@ -1,6 +1,7 @@
 import numpy as np
 from numpy.testing import assert_allclose
 
+import magpylib as magpy
 from magpylib._src.fields.field_BH_cuboid import magnet_cuboid_field
 from magpylib._src.fields.field_BH_cylinder_segment import magnet_cylinder_segment_field
 from magpylib._src.fields.field_BH_dipole import dipole_field
@@ -355,3 +356,192 @@ def test_field_line_v4():
         ]
     )
     np.testing.assert_allclose(B, Btest)
+
+
+def test_triangle1():
+    """test core triangle VS cube"""
+    obs = np.array([(3, 4, 5)] * 4)
+    mag = np.array([(0, 0, 333)] * 4)
+    fac = np.array(
+        [
+            [(-1, -1, 1), (1, -1, 1), (-1, 1, 1)],  # top1
+            [(1, -1, -1), (-1, -1, -1), (-1, 1, -1)],  # bott1
+            [(1, -1, 1), (1, 1, 1), (-1, 1, 1)],  # top2
+            [(1, 1, -1), (1, -1, -1), (-1, 1, -1)],  # bott2
+        ]
+    )
+    b = magpy.core.triangle_field("B", obs, mag, fac)
+    b = np.sum(b, axis=0)
+
+    obs = np.array([(3, 4, 5)])
+    mag = np.array([(0, 0, 333)])
+    dim = np.array([(2, 2, 2)])
+    bb = magpy.core.magnet_cuboid_field("B", obs, mag, dim)[0]
+
+    np.testing.assert_allclose(b, bb)
+
+
+def test_triangle2():
+    """test core single triangle vs same surface split up into 4 triangles"""
+    obs = np.array([(3, 4, 5)])
+    mag = np.array([(111, 222, 333)])
+    fac = np.array(
+        [
+            [(0, 0, 0), (10, 0, 0), (0, 10, 0)],
+        ]
+    )
+    b = magpy.core.triangle_field("B", obs, mag, fac)
+    b = np.sum(b, axis=0)
+
+    obs = np.array([(3, 4, 5)] * 4)
+    mag = np.array([(111, 222, 333)] * 4)
+    fac = np.array(
+        [
+            [(0, 0, 0), (3, 0, 0), (0, 10, 0)],
+            [(3, 0, 0), (5, 0, 0), (0, 10, 0)],
+            [(5, 0, 0), (6, 0, 0), (0, 10, 0)],
+            [(6, 0, 0), (10, 0, 0), (0, 10, 0)],
+        ]
+    )
+    bb = magpy.core.triangle_field("B", obs, mag, fac)
+    bb = np.sum(bb, axis=0)
+
+    np.testing.assert_allclose(b, bb)
+
+
+def test_triangle3():
+    """test core tetrahedron vs cube"""
+    ver = np.array(
+        [
+            [(1, 1, -1), (1, 1, 1), (-1, 1, 1), (1, -1, 1)],
+            [(-1, -1, 1), (-1, 1, 1), (1, -1, 1), (1, -1, -1)],
+            [(-1, -1, -1), (-1, -1, 1), (-1, 1, -1), (1, -1, -1)],
+            [(-1, 1, -1), (1, -1, -1), (-1, -1, 1), (-1, 1, 1)],
+            [(1, -1, -1), (1, 1, -1), (1, -1, 1), (-1, 1, 1)],
+            [(-1, 1, -1), (-1, 1, 1), (1, 1, -1), (1, -1, -1)],
+        ]
+    )
+
+    mags = [
+        [1.03595366, 0.42840487, 0.10797529],
+        [0.33513152, 1.61629547, 0.15959791],
+        [0.29904441, 1.32185041, 1.81218046],
+        [0.82665456, 1.86827489, 1.67338911],
+        [0.97619806, 1.52323106, 1.63628455],
+        [1.70290645, 1.49610608, 0.13878711],
+        [1.49886747, 1.55633919, 1.41351862],
+        [0.9959534, 0.62059942, 1.28616663],
+        [0.60114354, 0.96120344, 0.32009221],
+        [0.83133901, 0.7925518, 0.64574592],
+    ]
+
+    obss = [
+        [0.82811352, 1.77818627, 0.19819379],
+        [0.84147235, 1.10200857, 1.51687527],
+        [0.30751474, 0.89773196, 0.56468564],
+        [1.87437889, 1.55908581, 1.10579983],
+        [0.64810548, 1.38123846, 1.90576802],
+        [0.48981034, 0.09376294, 0.53717129],
+        [1.42826412, 0.30246674, 0.57649909],
+        [1.58376758, 1.70420478, 0.22894022],
+        [0.26791832, 0.36839769, 0.67934335],
+        [1.15140149, 0.10549875, 0.98304184],
+    ]
+
+    for mag in mags:
+        for obs in obss:
+            obs6 = np.tile(obs, (6, 1))
+            mag6 = np.tile(mag, (6, 1))
+            b = magpy.core.magnet_tetrahedron_field("B", obs6, mag6, ver)
+            h = magpy.core.magnet_tetrahedron_field("H", obs6, mag6, ver)
+            b = np.sum(b, axis=0)
+            h = np.sum(h, axis=0)
+
+            obs1 = np.reshape(obs, (1, 3))
+            mag1 = np.reshape(mag, (1, 3))
+            dim = np.array([(2, 2, 2)])
+            bb = magpy.core.magnet_cuboid_field("B", obs1, mag1, dim)[0]
+            hh = magpy.core.magnet_cuboid_field("H", obs1, mag1, dim)[0]
+            np.testing.assert_allclose(b, bb)
+            np.testing.assert_allclose(h, hh)
+
+
+def test_triangle4():
+    """test core tetrahedron vs cube"""
+    obs = np.array([(3, 4, 5)] * 6)
+    mag = np.array([(111, 222, 333)] * 6)
+    ver = np.array(
+        [
+            [(1, 1, -1), (1, 1, 1), (-1, 1, 1), (1, -1, 1)],
+            [(-1, -1, 1), (-1, 1, 1), (1, -1, 1), (1, -1, -1)],
+            [(-1, -1, -1), (-1, -1, 1), (-1, 1, -1), (1, -1, -1)],
+            [(-1, 1, -1), (1, -1, -1), (-1, -1, 1), (-1, 1, 1)],
+            [(1, -1, -1), (1, 1, -1), (1, -1, 1), (-1, 1, 1)],
+            [(-1, 1, -1), (-1, 1, 1), (1, 1, -1), (1, -1, -1)],
+        ]
+    )
+    b = magpy.core.magnet_tetrahedron_field("B", obs, mag, ver)
+    b = np.sum(b, axis=0)
+
+    obs = np.array([(3, 4, 5)])
+    mag = np.array([(111, 222, 333)])
+    dim = np.array([(2, 2, 2)])
+    bb = magpy.core.magnet_cuboid_field("B", obs, mag, dim)[0]
+
+    np.testing.assert_allclose(b, bb)
+
+
+def test_triangle5():
+    """special case tests on edges - result is continuous and 0 for vertical component"""
+    btest1 = [
+        [26.29963526814195, 15.319834473660082, 0.0],
+        [54.91549594789228, 41.20535983076747, 0.0],
+        [32.25241487782939, 15.087161660417559, 0.0],
+        [10.110611199952707, -11.41176203622237, 0.0],
+        [-3.8084378251737285, -30.875600143560657, -0.0],
+        [-15.636505140623612, -50.00854548249858, -0.0],
+        [-27.928308992688645, -72.80800891847107, -0.0],
+        [-45.34417750711242, -109.5871836961927, -0.0],
+        [-36.33970306054345, 12.288824457077656, 0.0],
+        [-16.984738462958845, 4.804383318447626, 0.0],
+    ]
+
+    btest2 = [
+        [15.31983447366009, 26.299635268142033, 0.0],
+        [41.20535983076747, 54.91549594789104, 0.0],
+        [-72.61316618947018, 32.25241487782958, 0.0],
+        [-54.07597251255013, 10.110611199952693, 0.0],
+        [-44.104089712909634, -3.808437825173785, -0.0],
+        [-36.78005591314963, -15.636505140623605, -0.0],
+        [-30.143798442143236, -27.92830899268858, -0.0],
+        [-21.886855846306176, -45.34417750711366, -0.0],
+        [12.288824457077965, -36.33970306054315, 0.0],
+        [4.80438331844773, -16.98473846295874, 0.0],
+    ]
+
+    n = 10
+    ts = np.linspace(-1, 6, n)
+    obs1 = np.array([(t, 0, 0) for t in ts])
+    obs2 = np.array([(0, t, 0) for t in ts])
+    mag = np.array([(111, 222, 333)] * n)
+    ver = np.array([[(0, 0, 0), (0, 5, 0), (5, 0, 0)]] * n)
+
+    b1 = magpy.core.triangle_field("H", obs1, mag, ver)
+    np.testing.assert_allclose(btest1, b1)
+    b2 = magpy.core.triangle_field("H", obs2, mag, ver)
+    np.testing.assert_allclose(btest2, b2)
+
+
+def test_triangle6():
+    """special case tests on corners - result is nan"""
+    obs1 = np.array([(0, 0, 0)])
+    obs2 = np.array([(0, 5, 0)])
+    obs3 = np.array([(5, 0, 0)])
+    mag = np.array([(111, 222, 333)])
+    ver = np.array([[(0, 0, 0), (0, 5, 0), (5, 0, 0)]])
+    b1 = magpy.core.triangle_field("B", obs1, mag, ver)
+    b2 = magpy.core.triangle_field("B", obs2, mag, ver)
+    b3 = magpy.core.triangle_field("B", obs3, mag, ver)
+
+    for b in [b1, b2, b3]:
+        np.testing.assert_equal(b, np.array([[np.nan] * 3]))
diff --git a/tests/test_obj_Facet.py b/tests/test_obj_Facet.py
new file mode 100644
index 000000000..19784d5ea
--- /dev/null
+++ b/tests/test_obj_Facet.py
@@ -0,0 +1,86 @@
+import numpy as np
+
+import magpylib as magpy
+from magpylib._src.exceptions import MagpylibMissingInput
+
+
+def test_Triangle_repr():
+    """Triangle repr test"""
+    line = magpy.misc.Triangle()
+    assert line.__repr__()[:8] == "Triangle", "Triangle repr failed"
+
+
+def test_triangle_input1():
+    """test obj-oriented triangle vs cube"""
+    obs = (1, 2, 3)
+    mag = (0, 0, 333)
+    vert = np.array(
+        [
+            [(-1, -1, 1), (1, -1, 1), (-1, 1, 1)],  # top1
+            [(1, -1, -1), (-1, -1, -1), (-1, 1, -1)],  # bott1
+            [(1, -1, 1), (1, 1, 1), (-1, 1, 1)],  # top2
+            [(1, 1, -1), (1, -1, -1), (-1, 1, -1)],  # bott2
+        ]
+    )
+    coll = magpy.Collection()
+    for v in vert:
+        coll.add(magpy.misc.Triangle(mag, v))
+    cube = magpy.magnet.Cuboid(mag, (2, 2, 2))
+
+    b = coll.getB(obs)
+    bb = cube.getB(obs)
+
+    np.testing.assert_allclose(b, bb)
+
+
+def test_triangle_input3():
+    """test core triangle vs objOriented triangle"""
+
+    obs = np.array([(3, 4, 5)] * 4)
+    mag = np.array([(111, 222, 333)] * 4)
+    vert = np.array(
+        [
+            [(0, 0, 0), (3, 0, 0), (0, 10, 0)],
+            [(3, 0, 0), (5, 0, 0), (0, 10, 0)],
+            [(5, 0, 0), (6, 0, 0), (0, 10, 0)],
+            [(6, 0, 0), (10, 0, 0), (0, 10, 0)],
+        ]
+    )
+    b = magpy.core.triangle_field("B", obs, mag, vert)
+    b = np.sum(b, axis=0)
+
+    tri1 = magpy.misc.Triangle(mag[0], vertices=vert[0])
+    tri2 = magpy.misc.Triangle(mag[0], vertices=vert[1])
+    tri3 = magpy.misc.Triangle(mag[0], vertices=vert[2])
+    tri4 = magpy.misc.Triangle(mag[0], vertices=vert[3])
+
+    bb = magpy.getB([tri1, tri2, tri3, tri4], obs[0], sumup=True)
+
+    np.testing.assert_allclose(b, bb)
+
+
+def test_empty_object_initialization():
+    """empty object init and error msg"""
+
+    fac = magpy.misc.Triangle()
+
+    def call_show():
+        """dummy function call show"""
+        fac.show()
+
+    np.testing.assert_raises(MagpylibMissingInput, call_show)
+
+    def call_getB():
+        """dummy function call getB"""
+        fac.getB()
+
+    np.testing.assert_raises(MagpylibMissingInput, call_getB)
+
+
+def test_barycenter():
+    """call barycenter"""
+    mag = (0, 0, 333)
+    vert = ((-1, -1, 0), (1, -1, 0), (0, 2, 0))
+    face = magpy.misc.Triangle(mag, vert)
+    bary = np.array([0, 0, 0])
+    np.testing.assert_allclose(face.barycenter, bary)
diff --git a/tests/test_obj_Tetrahedron.py b/tests/test_obj_Tetrahedron.py
new file mode 100644
index 000000000..0cfa35797
--- /dev/null
+++ b/tests/test_obj_Tetrahedron.py
@@ -0,0 +1,69 @@
+from doctest import REPORT_CDIFF
+
+import numpy as np
+
+import magpylib as magpy
+from magpylib._src.exceptions import MagpylibBadUserInput
+from magpylib._src.exceptions import MagpylibMissingInput
+
+
+def test_Tetrahedron_repr():
+    """Tetrahedron repr test"""
+    tetra = magpy.magnet.Tetrahedron()
+    assert tetra.__repr__()[:11] == "Tetrahedron", "Tetrahedron repr failed"
+
+
+def test_tetra_input():
+    """test obj-oriented triangle vs cube"""
+    obs = (1, 2, 3)
+    mag = (111, 222, 333)
+    vert_list = [
+        [(1, 1, -1), (1, 1, 1), (-1, 1, 1), (1, -1, 1)],
+        [(-1, -1, 1), (-1, 1, 1), (1, -1, 1), (1, -1, -1)],
+        [(-1, -1, -1), (-1, -1, 1), (-1, 1, -1), (1, -1, -1)],
+        [(-1, 1, -1), (1, -1, -1), (-1, -1, 1), (-1, 1, 1)],
+        [(1, -1, -1), (1, 1, -1), (1, -1, 1), (-1, 1, 1)],
+        [(-1, 1, -1), (-1, 1, 1), (1, 1, -1), (1, -1, -1)],
+    ]
+
+    coll = magpy.Collection()
+    for v in vert_list:
+        coll.add(magpy.magnet.Tetrahedron(mag, v))
+
+    cube = magpy.magnet.Cuboid(mag, (2, 2, 2))
+
+    b = coll.getB(obs)
+    bb = cube.getB(obs)
+    np.testing.assert_allclose(b, bb)
+
+    h = coll.getH(obs)
+    hh = cube.getH(obs)
+    np.testing.assert_allclose(h, hh)
+
+
+def test_tetra_bad_inputs():
+    """test obj-oriented triangle vs cube"""
+
+    bad_inputs = [
+        1,
+        [[(1, 1, -1), (1, 1, 1), (-1, 1, 1), (1, -1, 1)]] * 2,
+        [(1, 1, -1), (1, 1, 1), (-1, 1, 1)],
+        "123",
+    ]
+
+    mag = (111, 222, 333)
+    for bad in bad_inputs:
+
+        def test_function():
+            """bad vert input"""
+            magpy.magnet.Tetrahedron(mag, bad)
+
+        np.testing.assert_raises(MagpylibBadUserInput, test_function)
+
+
+def test_tetra_barycenter():
+    """get barycenter"""
+    mag = (111, 222, 333)
+    vert = [(1, 1, -1), (1, 1, 1), (-1, 1, 1), (1, -1, 1)]
+    tetra = magpy.magnet.Tetrahedron(mag, vert)
+    np.testing.assert_allclose(tetra.barycenter, (0.5, 0.5, 0.5))