A Computational Approach to Magnetic Force Feedback Design
Masa Ogata

Yuki Koyama

National Institute of Advanced Industrial Science andTechnology (AIST)Japan

National Institute of Advanced Industrial Science andTechnology (AIST)
Japan

(a)

(b)

(c)

(d)

Figure 1: Computational design of magnetic force feedback in digital fabrication. (a) The user specifies the target 3D model(consisting of a movable part and a fixed part) and the desired haptic force feedback. (b) Our method computationally solves
an inverse problem to obtain an optimal arrangement of magnets that renders the desired force feedback. (c) The user caneasily fabricate the model and embed magnets into the holes. (d) The user obtains the desired haptic force feedback during
interacting with the fabricated object.

ABSTRACT
We present a computational approach to haptic design embedded ineveryday tangible interaction with digital fabrication. To generate
haptic feedback, the use of permanent magnets as the mechanism
potentially contributes to simpleness and robustness; however, it isdifficult to manually design how magnets should be embedded in
the objects. Our approach enables the inverse design of magnetic
force feedback; that is, we computationally solve an inverse problem to obtain an optimal arrangement of permanent magnets that
renders the user-specified haptic sensation. To solve the inverse
problem in a practical manner, we also present techniques on magnetic simulation and optimization. We demonstrate applications to
explore the design possibility of augmenting digital fabrication foreveryday use.

CCS CONCEPTS• Human-centered computing → Human computer interaction (HCI).

KEYWORDS
Haptics, Magnetic force feedback, Inverse design, Optimization,
Digital fabrication
ACM Reference Format:
Masa Ogata and Yuki Koyama. 2021. A Computational Approach to Magnetic
Force Feedback Design. In CHI Conference on Human Factors in Computing
CHI ’21, May 8–13, 2021, Yokohama, Japan© 2021 Copyright held by the owner/author(s). Publication rights licensed to ACM.
This is the author’s version of the work. It is posted here for your personal use. Not
for redistribution. The definitive Version of Record was published in CHI Conferenceon Human Factors in Computing Systems (CHI ’21), May 8–13, 2021, Yokohama, Japan,
https://doi.org/10.1145/3411764.3445631.

Systems (CHI ’21), May 8–13, 2021, Yokohama, Japan. ACM, New York, NY,
USA, 11 pages. https://doi.org/10.1145/3411764.3445631

1

INTRODUCTION

Haptic design is essential for humans to understand the functionality of tools. In tangible interaction, producing appropriate haptics
can enhance usability and a sense of immersion. Actuators cancontribute to haptic design with interactivity; however, the scope islimited considering the implementation constraint. Aiming designfor everyday use, we focus on a physical mechanism for haptic
design. Some of the mechanical components are designed as a haptic device, but its design is generally normalized in industry andnot flexible. The magnetic force can be a mechanism for rendering
haptics and has a significant advantage with digital fabrication
since it is simple, compact, contactless, and non-abrasion. Hence,
our motivation is to achieve haptic design using magnetic forces
embedded in tangible interaction.
However, designing haptics using magnetic forces is highly challenging. Magnetic force is non-intuitive for the kinesthetic senseof the human body. For example, the attraction force between two
magnets changes unexpectedly in close distance. Besides, for our
purpose, a magnetic force can be perceived when the user moves
the movable component of an everyday object with embedded magnets; in this case, the magnetic force is hugely affected by not only
the trajectory of the moving magnet but also the orientations ofother magnets embedded in the fixed component. With some otherdifficulties, the user will face problems when deciding the spatial
relationship of all magnets because the manual design of magnets’
position and orientation in three-dimensional models is a highly
complicated task.
In this paper, we propose a computational approach to magnetic
force feedback design. As shown in Figure 1, the user effortlessly

CHI ’21, May 8–13, 2021, Yokohama, Japan

2 RELATED WORKS
2.1 Haptic Design by Magnetic Force

(a)

(b)

0

(c)
Container 1(C)builds a mechanical component for tangible interaction by specifying the desired curve of the force feedback by drawing. Even
with the support of the forward simulation [18], it was difficult
for the user to design by manual. To achieve our goal, we formulate this inverse design problem from the perspective of numerical
optimization and propose an optimization method with several
non-trivial techniques. We demonstrate the effectiveness of thisapproach by creating various objects with desired haptics using
digital fabrication.

Masa Ogata and Yuki Koyama

1(C)Container 2

Container 1
Trajectory C

Container 2

Figure 2: Input from the user. It consists of (a) the desired
curve of haptic sensation, (b) the geometry and maximum
number of magnets, and (c) the containers and trajectory.

2.3

Computational Inverse Design

Magnetic force has not been widely discussed as a potential sourceof haptic feedback. To the best of our knowledge, Hayward [9]firstly reported the relationship between haptics and magnetic
force, with mentioning similarity to lateral force fields reported
in [26]. As the capability of magnetic force, haptic design with amagnetic force is a popular method in human-computer interaction(HCI). Designing tactile feedback using a rubber magnet with amagnetization pattern is done by [35, 36]. Mechamagnets [39] performed design exploration of a haptic mechanism with permanent
magnets. Magneto-haptics [18] also explored forward simulation
for designing magnetic force with permanent magnets. However,a relationship between magnetic force and haptics has not been
discussed in previous works. In this paper, we formally discuss thisrelationship to solve the inverse design to achieve the automatic
design of the magnetic force.
Active haptic devices, such as actuator and vibrator, had beenused for haptics with computation from an early stage [14]. We
here mention some examples using the active haptic device witha magnetic force. MaglevHaptics [10] performed a haptic device
designed with magnetic levitation. FingerFlux [34] is a technicalapproach to generate haptic sensation via magnetic flux from anarray of the electromagnet. FluxPaper [19] tried vibration via magnetized paper using array of electromagnets. Omni [13] is a haptic
feedback system in which sensing and actuation are enabled bydetecting the 3D position of a magnet embedded in a passive tool.

Computational design methods aim at finding optimal design parameters that maximize desired properties by computationally solving optimization problems. Researchers have taken this approach invarious domains [3, 12, 20, 25, 32]; for example, the target property
could be the balancing capability of 3D printed objects [3, 25] or theuser performance in graphical user interface [20]. Inverse design
is considered as a variant of computational design, where the user
directly specifies the property that the final outcome should have,and the optimization solver tries to reproduce the property. Thisapproach is useful for fabrication scenarios, and researchers haveinvestigated inverse design methods that reproduce user-specified
mechanical properties [27, 38], deformations [4, 29], colors [2, 6],
and shapes [8, 15] of fabricated objects. Our work investigates aninverse design approach to reproduce the user-specified haptics
using embedded permanent magnets.
Several researchers have investigated the inverse design of haptics in different scenarios. Fujinawa et al. [7] proposed a methodof designing hand-held objects that provide illusional haptic shapeperception. Piovarci et al. [24] proposed a method of optimizing 3D
printing materials to design perceptual elasticity in touch interaction. Researchers have also investigated methods of achieving theinverse design of tactile textural feedback for fingertips [31] andstyli for drawing [22, 23]. Our work is the first to investigate theinverse design of magnetic force feedback, which enables end-users
to design new interactions with everyday objects, and we propose
several techniques to achieve this concept.

2.2

3

Simulation in Electromagnetism

There are mainly two approaches to simulate electromagnetism; analytical and numerical. For solving the magnetic force between two
magnets, some analytical calculations were proposed in the earlystage. As a numerical solution, the finite element method (FEM) is
the most famous for solving many problems in physics. However,the computation cost of FEM in three-dimensional space is veryhigh. An approximate analytical approach to calculate magnetic
forces between two magnets was developed using Taylor expansion [37]. This approach was also utilized in the computer graphics
community [30]. Magneto-Haptics [18] also utilized this approachfor forward simulation of magnetic force to estimate haptic feedback for designing tangible interaction. We follow them and use
this approximate approach. Note that some advanced magnetic
simulation methods tailored for specialized scenarios [11, 16] havebeen proposed recently, but these features are not necessary for ourtarget scenarios. Also note that the concepts from magnetostatics
used in this paper are explained in details in [5].

INTERACTION OVERVIEW

We propose a computational approach to design magnetic force
feedback in three-dimensional models for digital fabrication. A
spatial arrangement of magnets that fulfills the desired magnetic
force feedback can be automatically explored and optimized insoftware. The user specifies the desired curve of force feedback bydrawing and then obtains a functional three-dimensional modelwith spaces hollowed out for permanent magnets. The software only
requires some instruction about the movable component, without
asking any knowledge of magnetism.

3.1

Input and Output

The user needs to specify the followings as input (see Figure 2).
• The desired curve of haptic sensation (annotated as H (t)).
• The geometry of containers.
• The trajectory that one of the containers moves with.
• The geometry of magnet for each container.

A Computational Approach to Magnetic Force Feedback Design

CHI ’21, May 8–13, 2021, Yokohama, Japan

Table 1: Geometry of the magnets used in the applications. Sizes represent Width-Depth-Height for cubes, andDiameter-Height for cylinders (in millimeters). All magnets
are magnetized in the height direction.

Slider (1)Slider (2)Slider (3)Dial (1)Dial (2)Case (1)Container 1

Container 2

Type

Size

Type

Size

Cube
Cube
Cube
CylinderCylinder
Cube

4×3×6
5×5×2
4×4×4
4×4
4×3
3 × 24 × 2

CylinderCube
CylinderCylinderCylinderCylinder

3×6
5×5×2
3×3
4×4
3×4
3×6

• The maximum number of magnets.
The software then outputs an optimal arrangement of magnets, i.e.,
how we should embed magnets into the containers. In this work,
the number of containers is assumed to be two for simplicity. As forthe magnet material, the Neodymium magnet (N35 grade, which isvery popular in the market) is used by default. The user can change
the magnetic material by setting the value of residual magnetism.

3.2

Example Applications

Figure 3 outlines example applications of physical interface components that demonstrate haptic design by magnetic force. Threetypes of components are demonstrated: slider, dial, and case. The
user prepared the 3D models, and the system automatically hollowed holes for embedding magnets after fabrication by 3D printing.
The used magnets are summarized in Table 1. These geometries
were empirically chosen for each applications after several trials
with considering the solution quality.
Slider. This component has a small movable block that moves
on the linear trajectory. Inputting different kinds of target curves,
we obtained three unique arrangements. Slider (1) only uses threemagnets to achieve a single cycle of a sawtooth wave. Slider (2)performs a non-trivial arrangement to achieve linear force feedback
only using three magnets; besides the magnetic force between two
magnets is known as non-linear to the distance. Slider (3) shows thatthe system managed sharp spots by arranging magnets effectively.
Dial. This component has a movable disk-shape part that moves
on a rotational trajectory. The rotation angle is limited to 180 degrees on the system. For Dial (1), we tried to find the optimal arrangements which achieve bump and hole illusion. Dial (2) performsa small spot to give a stable position of haptic sensation.
Case. This component has single rid moves on a rotational trajectory with two holes. Usually, the magnetic force implemented for
those case design is too abrupt for the user when closing. We tried
to create linear force feedback to ease the magnetic force when
closing the rid. In order to reduce the number of magnets, the flatand long magnet for container 1 was selected. The rotation angle
is limited to 45° on the system.

3.3

Difficulty of Problem

The need for inverse design comes from the difficulties in manual
design. There are mainly two reasons. (1) Small changes in theorientation of any magnet drastically affect the magnetic force.
This is because the magnetic force in the magnetic field is positiondependent, and the superposed magnetic field leaked from multiple
magnets is highly complicated. (2) The space of parameters to
explore is highly broad and indirect. The magnetic force applied tothe movable magnet is varying on the trajectory, which requires
careful consideration of how multiple magnets affect the timevarying magnetic force in combination.

4

TECHNICAL OVERVIEW

We formulate the haptics design problem as an inverse optimization
problem. That is, given a desired haptic potential field (describedas haptic sensation in Section 3) specified by the user, we wantto find an optimal arrangement of magnets that reproduces thedesired haptic potential field as closely as possible. The magnetic
simulation is iteratively performed in our optimization process.
Figure 4 shows the overview of the process.
Notations. We denote by H ∈ H a haptic potential field (which
can be represented as a curve), where H is the set of all the possible
haptic potential fields. We especially denote by H target ∈ H thetarget haptic potential field that is specified by the user. We denoteby x ∈ X an arrangement of magnets, where X is the set of all thepossible magnet arrangements. Note that the number of magnets
is not a constant in our problem setting (i.e., it can be an arbitrary
number that is more than or equal to two), and thus we use x as atrans-dimensional variable; for example, x may implicitly represent
a spatial arrangement of two magnets, or that of five magnets. We
write H (x) ∈ H for the simulated haptic potential field that themagnet arrangement x generates. Finally, we denote by x∗ ∈ X the
optimal magnet arrangement that we want to obtain by solving the
optimization problem.
Mathematical Formulation. This inverse design problem can bemathematically described as an optimization problem:x∗ = arg min d(H (x), H target )x∈X
subject toall magnets are inside the container,
all magnets are collision-free(1)

where d : H × H → R ≥0 1 is a distance function that quantifies
how much the two haptic potential fields differ (formally definedin Section 7.2). For example, d(H 1 , H 2 ) returns zero when H 1 andH 2 are exactly the same, and it returns a positive value when H 1
is different from H 2 . The solution x∗ is the magnet arrangementthat minimizes the distance between the simulated haptic potential
field and the target haptic potential field, while it satisfies the twoconstraints, “all magnets are inside the container” and “all magnets
are collision-free.”Technical Challenges. This optimization problem is highly challenging to solve; it is considered as a constrained non-linear global
optimization, and the number of its dimensions is not fixed since
1 Here R≥0 denotes the set of real numbers that are greater than or equal to zero. We
will use similar notations in the reminder of the paper.

CHI ’21, May 8–13, 2021, Yokohama, Japan

Masa Ogata and Yuki Koyama

3D ModelH(t)[N · m]

Found optimal arrangement

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

H(t)[N · m]

Target / reproduced curve

H(t)[N · m]

H(t)[N · m]

Physical result

Slider (1)Slider (2)Slider (3)Dial (1)Dial (2)Case (1)Figure 3: 3D models and result of optimization with actual output fabricated by 3D printer and embeded permanent magnets.
Target haptic
potential ﬁeld

Optimization solver (Sec. 7)
Magnetarrangement

Max num.
magnets

Optimal magnet
arrangement

Hapticpotential ﬁeld

Magnetic simulator (Sec. 6)

Figure 4: Overview of the inverse design process. The input to the optimization solver consists of the target haptic potential field H target and the maximum number of magnets n max . The optimization solver then performs iterative search, which
involves the magnetic simulator in its loop. Finally, the optimization solver outputs the optimal magnet arrangement x∗ .
we do not know beforehand how many magnets should be used. We
will detail the technical challenges in Section 7.1. These challengings prevent us from directly applying off-the-shelf optimization
libraries; they do not accept this problem type, or at best, they will
fail to find plausible solutions in practical amount of time. This
motivates us to develop techniques to make this problem tractable.
Technical Contributions. In addition to the primary contribution,i.e., the proposition and demonstration of the computational approach to magnetic force feedback design, we also offer the following technical contributions to achieve this approach. In Section 5,
we formalize the definition of the haptic potential field and also
describe how it can be interpreted from the viewpoint of physics.
In Section 6, we describe new techniques to efficiently simulatemagnetic forces to calculate the haptic potential field H (x). In Section 7, we describe how the above optimization problem can beapproximately solved by making non-trivial modifications to thisproblem.

5

DEFINITION OF HAPTIC POTENTIAL FIELD

This section defines and discusses the haptic potential field, denoted as H . The concept itself was firstly introduced by Ogata [18];however, the original paper lacks formal discussions in physics, andit only provides how the concept is implemented and how visually
intuitive it is. For example, although it has “potential” in its name,the paper does not discuss how it relates to the concept of potential
in physics. We describe the concept more formally and generally
from the viewpoint of physics so that haptics researchers who workon magnetism or want to extend the concept beyond magnetism
can refer to this section in the future.
As an intuition, we describe the haptic potential field (which
will be defined in Equation 8) in relation to the work that thefinger does while manipulating the movable container along withthe predetermined trajectory (as we will see in Equations 9 and
14). The finger perceives this work as force feedback during themanipulation.

A Computational Approach to Magnetic Force Feedback Design

5.1

CHI ’21, May 8–13, 2021, Yokohama, Japan

AssumptionsFirst, we assume that all motions occur in the quasi-static manner. That is, any time derivatives (e.g., velocity and acceleration)in the environment are considered sufficiently small and thus ignorable. This assumption also implies that the magnetic field in
the environment is always considered steady, and thus we canuse magnetostatics. We also assume that the movable container issufficiently small, and it is considered a point mass (with a massm ∈ R >0 ). Finally, we assume that the movable container can only
move along the specified trajectory C because of the design ofcontainers.

5.2

Trajectory

We parametrize the three-dimensional trajectory C by a parametert ∈ [a, b], where a represents the starting endpoint and b represents
the ending endpoint. Without loss of generality, we assume a = 0
and b = 1. Let r be the function that maps t and the corresponding
three-dimensional position on the trajectory C. For example, r(0) ∈R3 and r(1) ∈ R3 are the positions of the starting and ending
endpoints of the trajectory C, respectively. Note that the parametert does not have a physical interpretation, but we call it “time” justfor convenience when no confusion can arise.

5.3

ForcesLet Ftotal (t) ∈ R3 be the total force that the movable container
receives at the time t during the motion along the trajectory C.
Similarly, let Ffinger (t), Fmagnet (t), Fgravity (t), and Fcontainer (t) be
the force applied by the finger, the force applied by the magnetic
field, the force applied by the gravity, and that force applied by theother (fixed) container, respectively. Then, the relation,
Ftotal (t) = Ffinger (t) + Fmagnet (t) + Fgravity (t) + Fcontainer (t), (2)holds at any time t. By applying Newton’s second law to the movable container, we have Ftotal (t) = ma(t), where a(t) ∈ R3 is theacceleration of the movable container at the time t. By the quasistatic assumption, we have a(t) = 0 at any time t. Therefore,
Ffinger (t) = −Fmagnet (t) − Fgravity (t) − Fcontainer (t)(3)

holds at any time t.

5.4

Work Done by Finger

We now consider the work done by the application of the force tothe movable container by the finger. We denote by Wfinger (t) ∈ R
the work done until reaching the position r(t) from the starting
position r(0). By the definition of work, it is written as∫ t
Wfinger (t) =Ffinger (t ′ ) · r ′ (t ′ )dt ′,(4)0∂ r(t), and the centered dot represents the inner
where r ′ (t) = ∂tproduct between vectors.

5.5

Potential of Magnetic Force

The magnetic force, Fmagnet , needs a special care. Suppose an ideal
case of a uniformly magnetized magnet; it can be treated as a simple
magnetic moment. However, as actual magnets have various shapes

and magnetic fields created by other magnets are not uniform, magnets cannot be treated as magnetic moments to simulate Fmagnet
in a feasible way. Thus, we need to utilize numerical analysis tocalculate magnetic fields and the magnetic force.
Nevertheless, the force Fmagnet can be considered a conservative
force since it is fully specified by the spatial position of the movable container by the magnetostatics assumption. Thus, though nosimple closed form is available, we can still define a potential fieldUmagnet (t) such that Fmagnet (t) = −∇Umagnet (t).

5.6

Haptic Potential Field

Considering Equation 3 and Equation 4, the work Wfinger (t) can be
virtually decomposed into the component by conservative forcescon. (t), and the com(i.e., −Fmagnet (t) − Fgravity (t)), denoted by Wfinger
ponent by non-conservative forces (i.e., −Fcontainer (t)), denoted by
non-con. (t). That is,
Wfingercon.
non-con.
Wfinger (t) = Wfinger(t) + Wfinger
(t).

(5)

Let Wmagnet (t) be the work done by the application of the magnetic
force Fmagnet (t) to the movable container. Similarly, let Wgravity (t)be the work done by the application of the gravitational force
Fgravity (t) to the movable container. Using these notations, we havecon.
Wfinger(t) = −Wmagnet (t) − Wgravity (t).

(6)

Since they are the work done by conservative forces, it can be
written using their potentials as
con.
Wfinger(t) = ∆0→t Umagnet + ∆0→t Ugravity ,(7)

where ∆a→b U = U (b) − U (a). Now, we define the haptic potential
field2 H : [0, 1] → R as
H (t) = Umagnet (t) + Ugravity (t).

(8)

Since both components are potential fields, H is also a potential
field. Using this definition, Equation 7 is written ascon.
Wfinger(t) = ∆0→t H = H (t) − H (0),(9)

where H (0) is an arbitrary constant offset. Thus, the landscape (i.e.,the curve) of the haptic potential field can be interpreted as the plotof the work done by the finger against conservative forces (with anoffset).
It is noteworthy that the haptic potential field does not involve
any path-dependent line integrals in its definition; the landscapeof H is thus independent of the parameterization of the trajectory
C. This is why we call it a field. In other words, it is beneficial forthe definition to be independent of the force applied by the fixedcontainer, Fcontainer , since it is a non-conservative force (unlike
the gravity force) and thus prohibits us from modeling the haptic
sensation from the viewpoint of the potential field. Also, the force
may depend on many factors such as the choice of materials, which
requires additional assumptions to define the interaction scenario.
2 By abuse of notation, we here use H as the function of t (representing the value at apoint in the field). As described in Section 4, we also continue to use H as the functionof x (representing the entire field) when no confusion can arise.

CHI ’21, May 8–13, 2021, Yokohama, Japan

5.7

Masa Ogata and Yuki Koyama

Discussion: The Effect of Gravity

If the trajectory C involves only horizontal motions, then Ugravity (t)is constant at any time t, and thus ∆Ugravity = 0 holds. In such cases,
we can just ignore the effect of gravity in terms of the haptic potential field. Even when the trajectory C involves some up-down
motions, we can consider ∆Ugravity = 0 if the gravitational force
Fgravity = mg, where g ∈ R3 is the gravitational acceleration, issufficiently smaller than the magnetic force Fmagnet . In the following sections of this paper, we assume that we can ignore thegravitational effect for simplicity.

5.8

Discussion: The Effect of Friction

The force applied by the fixed container, Fcontainer (t), can be decomposed into the frictional component Ffriction (t), which aligns
with the moving direction, and the normal component Fnormal (t),which is orthogonal to the frictional component. That is,
Fcontainer (t) = Ffriction (t) + Fnormal (t).

(10)

Using this notation, the non-conservative component of the work,
non-con. , can be written as
Wfinger∫ t
non-con.
Wfinger(t) =(−Fcontainer (t ′ )) · r ′ (t ′ )dt ′
(11)
0
∫ t

=−
Ffriction (t ′ ) · r ′ (t ′ )
0

=−

∫ t
0

+ Fnormal (t ′ ) · r ′ (t ′ ) dt ′(12)

Ffriction (t ′ ) · r ′ (t ′ )dt ′,(13)

since Fnormal (t ′ ) · r ′ (t ′ ) is zero at any time t ′ .
If the frictional force Ffriction is sufficiently small (e.g., by choosnon-con. =ing smooth materials for fabrication), we can consider Wfinger
0. Thus, we have
Wfinger (t) = ∆0→t H = H (t) − H (0).

(14)

This relationship provides an additional interpretation on the haptic
potential field; the landscape (i.e., the curve) of the haptic potential
field can be interpreted as the plot of the work done by the finger(with an offset) in an frictionless environment.

6

MAGNETIC SIMULATION TECHNIQUESWe describe both the numerical analysis method to calculate themagnetic force F applied to the magnet moved by the finger, and theformulation to calculate the haptic potential field H (t). We utilize
the Taylor expansion approach, which approximately calculates themagnetic force between two magnets. This approach was originallyproposed in [37] and used for simulation in computer graphics [30].
In the Taylor expansion approach, the mesh resolution contributes to the precision. While reducing the number of meshingblocks contributes to the efficiency of calculation cost, a naïve way
would sacrifice the precision. Thus, we also describe some techniques to improve the calculation cost while maintaining precision.

6.1

Magnetic Force by Taylor Expansion

Basic formula of magnetic force applied to a magnetic pole is written
as F = ∇(m · B). Here, B means magnetic flux density, which is

Figure 5: Magnetic force simulation with dipoles. The force
F that the magnet M A receives from the magnet M B (Left)is considered as the summation of forces between dipoles(Right).
strength of magnetic field though certain surface depending onmaterial. In our purpose, B is determined as leaked magnetic flux
from other magnets. To approximate the magnetic force calculation,
we apply Taylor expansion to small sections of the magnet, treating
as small magnetic moments (magnetic dipoles) as illustrated in
Figure 5. Let F be the force that the magnet M A receives from themagnet M B , and n A and n B be the number of grid mesh cells (i.e.,
dipoles) of M A and M B , respectively. Then, the force is calculated
as
nA
nA Õ
nB
Õ
Õ
F=
Fk =
Fi,k ,(15)
k =1

k =1 i=1

where Fk is the force that the k-th dipole of the magnet M A receives
from the entire magnet M B , and Fi,k is the force that the k-th dipoleof the magnet M A receives from the i-th dipole of the magnet M B .
Figure 5 illustrates this concept. The force Fi,k is calculated ashµ0
1·−15ni,k (mk · ni,k )(mi · ni,k )Fi,k =
44π ∥rk − ri ∥i
+ 3n(mk · mi ) + 3 mk (mi · ni,k ) + mi (mk · ni,k )(16)
where ni,k = (rk − ri )/∥rk − ri ∥, r is a position of magnetic dipole,and µ 0 is the magnetic constant.
As the magnetic force is superpositional, the magnetic force frommultiple magnets can be solved by summarizing each F with eachmagnet. Still, there is a problem of determining eligible n A and n B ,which is a split number of the magnet into dipoles. These numbers
significantly affect the computation cost. We describe a techniqueto adaptively determine the split number in the next subsection.

6.2

Haptic Potential Field

As discussed in Section 5.6, we assume ∆H = ∆U with ignoringgravitational force. Therefore, the haptic potential field H (t) is
treated as equal to the potential field of magnetic force U (t). The
potential field U (t) can be obtained from the formula of ∆U :∫ t
∆U = −
F(r(t ′ )) · r ′ (t ′ )dt ′ .
(17)0Note that, to extract the component of the magnetic force applied
to the direction along the trajectory, the magnetic force must be“normalized” with the direction vector r ′ (t) by inner product, as
indicated by this equation. By assuming that the trajectory is discretized into some sections, we calculate F(r(t)), r ′ (t), and their

A Computational Approach to Magnetic Force Feedback Design

CHI ’21, May 8–13, 2021, Yokohama, Japan

Algorithm 1: The outer-loop part of the inverse design.

Figure 6: Adaptive meshing. From left to right, the upper
magnet comes closer to the bottom magnet. The resolutionof the meshes adaptively changes based on the distance, bywhich the computational cost of the magnetic force simulation can be reduced while the accuracy remains sufficient.
inner product at each section of t. Then, the potential field U (t)can be calculated as a cumulative product3 of normalized magnetic
forces and divided lengths of these sections.

6.3

Adaptive Meshing

When calculating the magnetic force, a pair of closer dipoles affects
more than distant ones. To increase the precision of the simulation,
it is desirable to execute meshing the closer pair of dipoles finely.
However, fine meshing causes high computation costs. Adaptive
meshing, i.e., selecting blocks to separate into finer blocks, will
help the problem of computation costs. Following the previous
work [37], we decided on the condition to apply adaptive meshing.
Condition is determined by ρ ≥ 1/7, where ρ = d/r , r is distance
between two dipoles, and d is length of dipole. Figure 6 shows timeseries of adaptive meshing up to depth 4, which is the maximum
depth level in our implementation. The number of depth may affectthe precision of the magnetic force.

6.4

Input: The target haptic potential field, H target , and the
maximum number of magnets, n max
Output: An optimal magnet arrangement (i.e., x∗ )target ;1 Perform pre-processing to H
2 Insert the first magnet (i.e., M 1 ) to the movable container;
max do3 for i = 2, 3, . . . , n
4
Search for a good initial location and orientation to add
a new magnet to the fixed container;
5
Insert a new magnet (i.e., Mi ) to the found initial
location and orientation;
6
Perform the inner-loop procedure (Algorithm 2) toadjust the arrangement of the magnets {M 1 , . . . , Mi };7
if Termination condition is met then break;
8 end
9 return The arrangement of the magnets {M 1 , . . . , M i };Interpolation

To reduce computation cost, we also operate interpolation on themagnetic force. We multiplied section four times from the original
curve of magnetic force. For example, if the original section of t is
20, it becomes 80, meaning that the computation cost becomes fourtimes smaller. Due to its characteristics, the shape of the magnetic
force is always smooth. Therefore, applying the spline interpolationis reasonable, and it does not cause much deviation from the original
curve of the magnetic force.

7 OPTIMIZATION TECHNIQUES
7.1 Overview
We want to solve the optimization problem, formulated as Equation 1, in a practical manner. Even with the efficient magnetic force
simulation described in Section 6, this problem is still very challenging to solve. This is because the problem has the following
properties.
• Trans-dimensional: the number of magnets (and thus the
optimization dimensionality) is not fixed beforehand andneeds to be determined by the solver.
• Constrained: the optimization variables need to not only
minimize the objective but also satisfy the constraints.
3 https://numpy.org/doc/stable/reference/generated/numpy.cumprod.html• High-dimensional: the number of optimization variables islarge.
• Global: there can be many local minima, most of which are
not useful.
• Non-linear: both the objective and constraints are not simple
linear functions.
• Non-differentiable: the derivatives of the objective and constraint functions are not available.
To our knowledge, no off-the-shelf solvers can handle this problem.
We propose the following three approximation techniques tosolve the problem practically.
Approximation 1: Incrementally insert magnets one by one
and stop this procedure when a termination condition is met,
rather than simultaneously optimizing the arrangement ofmagnets and the number of magnets. This approximation
resolves the issue of the optimization problem being transdimensional.
Approximation 2: Sequentially solve a small (low-dimensional)
optimization problem for each magnet, rather than simultaneously solving a large (high-dimensional) optimization
problem for all magnets. This approximation resolves theissue of the optimization problem being high-dimensional.
Approximation 3: Apply the penalty method [17] to reformulate the constrained optimization problem into a series of
unconstrained optimization problems. This approximation
allows the use of unconstrained optimization algorithms.
Taking these three techniques into consideration, we describe theoverall optimization procedure as an “outer”-loop part, shown in
Algorithm 1, and an “inner”-loop part, shown in Algorithm 2. The
outer-loop part involves Approximation 1, and it iteratively calls
the inner-loop part (line 6) in its loop. The inner-loop part involves
Approximations 2 and 3, and it iteratively solves a much easier
optimization problem (line 4), which is still non-linear and nondifferentiable but can be handled by off-the-shelf solvers if wechoose them appropriately.
In addition to these approximation techniques, we also propose
a technique to estimate a reasonable initial solution when inserting

CHI ’21, May 8–13, 2021, Yokohama, Japan

Masa Ogata and Yuki Koyama

Using these functions, Equation 1 can be rewritten as

Algorithm 2: The inner-loop part of the inverse design.

x∗ = arg min d(H (x), H target )Input: The target haptic potential field, H target , and theinitial arrangement of the i magnets, {M 1 , . . . , Mi }Output: An optimal magnet arrangement of the i magnets
penalty do
1 for j = 1, 2, . . . , n
2Solve the optimization problem to adjust thearrangement of the magnet Mi ;3for k = 1, 2, . . . , i do
4
Solve the optimization problem to adjust thearrangement of the magnet Mk ;5end
6
Increase the penalty coefficient σ ;7 end
8 return The arrangement of the magnets {M 1 , . . . , M i };x∈Xsubject to

7.4

Objective Formulation: Distance Metric for
Haptic Potential Fields

The objective of the optimization problem, Equation 1, is definedas the distance between the simulated haptic potential field andthe user-specified haptic potential field. As the definition of the
distance function d, we use the line integral of their distance ateach position. That is, given two haptic potential fields, H 1 and H 2 ,defined along the trajectory C, the distance is calculated by∫ 1∫∥H 1 (t) − H 2 (t)∥ ds =∥H 1 (t) − H 2 (t)∥ dt
d(H 1 , H 2 ) =
C

0

 
 
N
1 Õ
i
i
H1
,
≈
− H2
N + 1 i=0
N
N

(18)

where N is the number of samples to approximately calculate theline integral (N = 20 in our implementation).

7.3

Constraint Formulation: Validity of Magnet
Arrangements

As we described in Section 4, we consider two constraints thatshould be satisfied:• all magnets are inside of the containers, and• all magnets are collision-free.

(19)

More specifically, when the magnet arrangement is collision-free(i.e., satisfying the constraint), c collision returns a negative value that
represents the shortest distance among the magnets in the scene.
When there are collisions (i.e., violating the constraint), c collision
returns a positive value that represents the largest penetration
depth among the collided magnets in the scene. We design c container
similarly. Our implementation uses FCL [21] to calculate thesefunction values.

a new magnet (line 4), which helps find a good local optimum
without sticking to undesirable local optima. This initial solution
estimation allows the use of efficient local optimization algorithms
instead of prohibitively heavy global optimization algorithms (suchas Bayesian optimization [28]).
In the following subsections, we first describe more detailed formulations of the objective (Section 7.2) and constraints (Section 7.3)
in Equation 1. We then describe details of the techniques and steps
in Algorithms 1 and 2.

7.2

 containerc(x) ≤ 0.
c collision (x) ≤ 0

We introduce constraint functions, c container : X → R and c collision :X → R, to mathematically model these two constraints. We design
these functions such that they return negative values when themagnet arrangement satisfies the constraints, and they return positive values when the magnet arrangement violates the constraints.

Pre-Processing of User Input

Before starting to solve the optimization problem, we perform preprocessing to the user-specified haptic potential field H target , as inline 1 in Algorithm 1. This is because the user input may include
undesirable small noises due to mouse manipulation accuracy. Thus,
we first apply the Savitzky–Golay filter4 to H target .

7.5

Incremental Insertion of Magnets andTermination Condition

At the beginning of the inverse design procedure, both the movableand fixed containers do not contain any magnets. After the preprocessing, we insert the first magnet (i.e., M 1 ) to the movable
container (line 2 in Algorithm 1). Then, we begin to perform theouter-loop procedure. In each iteration, we insert a new magnet (i.e.,
Mi if it is the i-th magnet in the scene) to the fixed container (line 5
in Algorithm 1) and then perform the inner-loop procedure (line 6
in Algorithm 1) to adjust the arrangement of the so-far inserted
magnets. We stop the outer-loop procedure when a termination
condition is met (line 7 in Algorithm 1). For this condition, we
simply set a threshold for the acceptable distance, d threshold , andcheck whether the current distance between the simulated haptic
potential field and the user-specified one is smaller than d threshold ,or not. If it is true, then the outer-loop part terminates and returns
the solution obtained so far to the user (line 9 in Algorithm 1). Our
implementation adaptively sets d threshold based on H target ; more
specifically, the 5% of the value, maxt H target (t) − mint H target (t).

7.6

Estimation of Good Initial Solution

It is crucial to have a good initial solution in our problem setting
since our problem is a global optimization with many useless local
minima. To obtain a good initial arrangement of the first magnet,
M 1 , we perform the following procedure.
• Discretize the volumes of the movable and fixed containers
using uniform grids.
• For each grid cell in the movable container, calculate the line
integral of the distance to the closest grid cell in the fixedcontainer along the trajectory.
• Choose the grid cell in the movable container that has thesmallest value and call it G movable .
4 https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html

A Computational Approach to Magnetic Force Feedback Design• Use the center of G movable as the initial position.
• Use “N-up” (i.e., the N pole faces up) as the initial orientation.

To obtain a good initial arrangement of the i-th magnet, Mi (i ≥ 2),
we perform the following procedure.

• Find the parameter t that produces the maximum residual(i.e., the parameter t in which the difference between thetarget and current simulated haptic potential values become
the largest) and call it t ′ .
• Move the movable container to t = t ′ .
• Choose the grid cell that is closest to G movable from the fixedcontainer.
• Use the center of the chosen grid cell as the initial position.
• Use either “N-up” or “S-up” (determined based on the signof the residual) as the initial orientation.

7.7

Penalty Method

The penalty method [17] is one of the methods that solve a constrained optimization problem by decomposing it into a series of
unconstrained optimization subproblems. It iteratively refines thesolution by using the solution in the previous subproblem as theinitial solution in the next subproblem. It adds penalty terms to theoriginal objective, each of which represents how much the corresponding constraint is violated. That is, the value of the penalty
term is zero if the constraint is satisfied, and the value becomeslarger as the constraint gets violated. It increases the weight forthe penalty term each time (line 6 in Algorithm 2), so that the final
solution (mostly) satisfies the constraints. The j-th subproblem tosolve Equation 19 is described as
1x∗ = arg min d(H (x), H target )
h
x∈X
n
o
+ σ j д(c container (x)) + wд(c collision (x)) , (20)where д(x) = max(0, x)2 , h = maxt H target (t) − mint H target (t),σ j ∈ R >0 is the penalty weight at the j-th iteration of the penalty
method, and w is a constant parameter for balancing the constraints(empirically set to 1.0 × 10−3 ). Note that normalizing the distanceby h1 here is only optional but this can make the setting of thepenalty weights less dependent on the user input H target . Our implementation solves this subproblem twice (i.e., n penalty = 2 inAlgorithm 2) using σ1 = 5.0 × 104 and σ2 = 2.0 × 105 by default.

7.8

Magnet-by-Magnet Optimization

Each magnet has the degrees of freedom of a rigid body transformation, SE(3). For this, we use the Cartesian coordinates representation (which is three-dimensional) to indicate the translational component and the rotation vector representation5 (whichis three-dimensional) to indicate the rotational component. Thus,
the degrees of freedom of the k-th magnet, xk , is a 6-dimensional
vector in our implementation. Note that the reason for using therotation vector representation as optimization variables is that thisrepresentation is smooth and unbounded, is free from gimbal lock,

CHI ’21, May 8–13, 2021, Yokohama, Japanand does not require normalization or orthogonalization (unlike
the Euler angles, quaternion, and matrix representations).
As described in Approximation 2 in Section 7.1, we solve
Equation 20 for each magnet sequentially, not simultaneously. By
this technique, we avoid directly solving the 6i-dimensional optimization problem, where i is the number of magnets, and insteadsolve 6-dimensional optimization subproblems (see lines 3–5 in
Algorithm 2). Now that each optimization problem is sufficiently
tractable, we apply the Nelder–Mead method6 available in SciPy[33], which is a popular choice for non-linear, non-differentiable,
and reasonably local and low-dimensional optimization problems.
Note that the reason for solving the 6-dimensional optimization
problem for the new magnet first (line 2) is that the new magnet
is likely to be in a relatively unreasonable configuration comparedto other magnets, and thus it is effective to adjust its arrangement
first.

8

TECHNICAL EVALUATION

Figure 3 shows the results of our inverse design procedure withvarious inputs, including haptic curves with linear slopes or multiple peaks, and trajectories with a straight or curved path. Evenfor these challenging inputs, our solver could successfully find reasonable and non-trivial magnet arrangements that well reproduce
the target haptic curves. Refer to the video figure as well, which
provides detailed visualizations of how our optimization proceduresearched for the magnet arrangements for these cases. These results
successfully validate that our proposed computational approach tomagnetic force feedback design can be achieved by our proposed
techniques.
We also conducted additional experiments to further validate
the effectiveness of each technique. Our implementation runs on
a MacBook Pro (Processor 2.6 GHz Quad-Core Intel Core i7 andMemory 16 GB). It is written with Python for optimization andC++ for magnetic simulation. The simulation part is accelerated byusing OpenCL and runs on GPU (AMD Radeon Pro 450).
Effectiveness of Initial Solution Estimation. Our optimization solver
has a heuristic initial solution estimation (Section 7.6) in its procedure. To validate the effectiveness of this estimation, we performeda custom optimization procedure that is identical to the proposed
one except that it randomly initializes the newly added magnet’s
position. For this, we used the same input as Dial (1) in Figure 3.
Figure 7 shows the result, in which the solver clearly failed to finda reasonable solution, and the reproduced haptic curve was quite
different from the target one. This implies that our target optimization problem has many useless local minima, and our technique of
initial solution estimation could reasonably resolve this issue.
Effectiveness of Magnet-by-Magnet Optimization. As discussed
in Section 7.8, we avoid directly solving a 6n-dimensional problem,
where n is the number of magnets, but we solve 6-dimensional
subproblems iteratively. To validate the effectiveness of this technique, we compared our optimization procedure with a custom
one in which the sequential 6-dimensional subproblems were replaced with a single 6n-dimensional problem. For this, we used
the same input as the second example in Figure 3. As a result, we

5 https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.transform.

Rotation.html

6 https://docs.scipy.org/doc/scipy/reference/optimize.minimize-neldermead.html

CHI ’21, May 8–13, 2021, Yokohama, Japan

Figure 7: Result of using random initialization instead of using our initialization technique. It failed to find a reasonable
solution.
observed that both the procedures could find reasonable solutions(more specifically, the magnet arrangements were different, but
the reproduced haptic curves were mostly identical). However, theprocedure that solved 6n-dimensional problems took 35 minutes
and 57 seconds to achieve convergence, while ours took only 3
minutes and 57 seconds.
Necessity of Penalty Method. We formulated the inverse design
problem as a constrained optimization problem, and our solver
handles the constraints by the penalty method (Section 7.7). To confirm the necessity of this technique, we performed an optimization
procedure without adding the penalty of the constraints. For this,we used the same input as Dial (1) in Figure 3. Although the target
haptic curve is well reproduced, two of the three magnets protrude
from the containers.

9

DISCUSSION

Physical and Perceptual Validation. We did not perform measurement of magnetic forces (e.g., using a force gauge) for validation ofmagnetic simulation. This is because Yuang et al. [37] already performed a comparison test with analytical calculation and validated
the simulation accuracy, and Thomaszewski et al. [30] also compared this approach [37] with the analytical approach [1]. As forthe perception, Ogata [18] already validated the perceptual aspectof haptic potential field.
Physically Invalid Inputs. It is possible that the user-specified
haptic potential field is physically invalid (i.e., not achievable byany magnet arrangements), such as the one with spikes or highfrequency waves. In such cases, the inverse design process cannot
converge to a satisfactory solution that well reproduces the target.
When this occurs, it is possible to warn the user that the target
haptics could not be achieved well and needs to be fixed.
Non Quasi-Static Settings. Our method assumes a quasi-static
setting as discussed in Section 5.1. This assumption is not valid
when, for example, the user manipulates the movable container
back and forth very quickly. In this case, the kinetic energy needs tobe considered in addition to the potential energies. It is an important
future work to model both the potential and kinetic energies as a
unified haptic model.
Frictions. The haptic potential field can be understood as thework done by the finger in a frictionless environment (Section 5.8).
This implies a limitation: our current formulation ignores the effectof frictions. We observed that the frictional force in our applications(Section 3.2) was not large compared to the magnetic force, and it did

Masa Ogata and Yuki Koyama

not qualitatively change the characteristics of the perceived haptic
curves. Thus, this limitation is not critical, at least in our cases.
Nonetheless, allowing users to design haptics in a friction-awaremanner would be an important future work, making the design
space broader. In terms of techniques, our theory (Section 5) isalready capable of quantifying the contribution of friction to haptic
feedback in a unified framework, and this goal can be achieved byintegrating the work done against the frictional force (Equation 13)
into the objective of the inverse design. Note that the frictional
force could be simulated by putting additional assumptions on thefriction model (e.g., Coulomb’s friction law), the direction of theforce applied by the fingers, and so on. Another challenge regarding
this future work is the user specification; the input haptic curve
may be less intuitive for users since it is no longer understood as afield (in a physics sense) with friction.
Gravity. Our current implementation focuses on the magnetic
potential as the source of force feedback. However, the definitionof the haptic potential field (Equation 8) allows us to design haptics
by taking gravitational potential into consideration as well. If we
handle trajectories that involve vertical movements, extending theimplementation using the gravitational effect would be useful tomodel the perceived haptics more accurately.
Optimization Solver Accuracy. We demonstrated that our optimization solver could produce useful and non-trivial magnet arrangements for various inputs. Yet, since it relies on approximation
techniques, it is possible that there exists a better solution than the
one found by the solver. For example, our solver cannot guaranteethat the number of used magnets is minimal.
Electromagnet. We focused on the use of permanent magnets(not electromagnets) since they need no power supply or electric
wiring and thus are suitable for embedded haptics in everyday objects. Nonetheless, it is an interesting future work to investigate
the use of electromagnets to render haptic potential fields. As electromagnets can be programmable, this would broaden the design
space of interactions. Note that our discussions on haptic potential
field (Section 5) is directly useful even in this advanced scenario.

10

CONCLUSION

We demonstrated that our approach could produce non-trivial designs to achieve desired magnetic force feedback, offering a new
way to augment physical components for effective tangible interaction. As the magnetic force is an essential part of the source of any
mechanical power, the combination with a computational approachwith inverse design has various potentials. We envision a futureworld where we design our living environment with the power ofcomputational design and physical simulation. With the current
development of digital fabrication technology, we believe the combination of computational design and haptic design backed up byphysics will open up new methodologies for interaction.

REFERENCES[1] G. Akoun and J. . Yonnet. 1984. 3D analytical calculation of the forces exerted
between two cuboidal magnets. IEEE Transactions on Magnetics 20, 5 (1984),
1962–1964. https://doi.org/10.1109/TMAG.1984.1063554

A Computational Approach to Magnetic Force Feedback Design[2] Thomas Auzinger, Wolfgang Heidrich, and Bernd Bickel. 2018. Computational
Design of Nanostructural Color for Additive Manufacturing. ACM Trans. Graph.
37, 4, Article 159 (July 2018), 16 pages. https://doi.org/10.1145/3197517.3201376[3] Moritz Bächer, Emily Whiting, Bernd Bickel, and Olga Sorkine-Hornung. 2014.
Spin-It: Optimizing Moment of Inertia for Spinnable Objects. ACM Trans. Graph.
33, 4, Article 96 (July 2014), 10 pages. https://doi.org/10.1145/2601097.2601157[4] Bernd Bickel, Moritz Bächer, Miguel A. Otaduy, Hyunho Richard Lee, Hanspeter
Pfister, Markus Gross, and Wojciech Matusik. 2010. Design and Fabrication ofMaterials with Desired Deformation Behavior. ACM Trans. Graph. 29, 4, Article
63 (July 2010), 10 pages. https://doi.org/10.1145/1778765.1778800[5] John MD Coey. 2010. Magnetism and magnetic materials. Cambridge university
press.
[6] Oskar Elek, Denis Sumin, Ran Zhang, Tim Weyrich, Karol Myszkowski, Bernd
Bickel, Alexander Wilkie, and Jaroslav Křivánek. 2017. Scattering-Aware Texture
Reproduction for 3D Printing. ACM Trans. Graph. 36, 6, Article 241 (Nov. 2017),
15 pages. https://doi.org/10.1145/3130800.3130890[7] Eisuke Fujinawa, Shigeo Yoshida, Yuki Koyama, Takuji Narumi, Tomohiro
Tanikawa, and Michitaka Hirose. 2017. Computational Design of Hand-Held
VR Controllers Using Haptic Shape Illusion. In Proceedings of the 23rd ACM
Symposium on Virtual Reality Software and Technology (Gothenburg, Sweden)(VRST ’17). Association for Computing Machinery, New York, NY, USA, Article
28, 10 pages. https://doi.org/10.1145/3139131.3139160[8] Ruslan Guseinov, Eder Miguel, and Bernd Bickel. 2017. CurveUps: Shaping
Objects from Flat Plates with Tension-Actuated Curvature. ACM Trans. Graph.
36, 4, Article 64 (July 2017), 12 pages. https://doi.org/10.1145/3072959.3073709[9] Vincent Hayward. 2008. A brief taxonomy of tactile illusions and demonstrations
that can be done in a hardware store. Brain research bulletin 75, 6 (2008), 742–752.
[10] Ralph Hollis. 2008. Maglev Haptics: Butterfly Haptic’s New User Interface
Technology. In ACM SIGGRAPH 2008 New Tech Demos (Los Angeles, California) (SIGGRAPH ’08). ACM, New York, NY, USA, Article 23, 1 pages. https://doi.org/10.1145/1401615.1401638[11] Seung-Wook Kim, Sun Young Park, and Junghyun Han. 2018. Magnetization
Dynamics for Magnetic Object Interactions. ACM Trans. Graph. 37, 4, Article 121(July 2018), 13 pages. https://doi.org/10.1145/3197517.3201402[12] Yuki Koyama and Masataka Goto. 2018. OptiMo: Optimization-Guided Motion
Editing for Keyframe Character Animation. In Proceedings of the 2018 CHI Conference on Human Factors in Computing Systems. Association for Computing
Machinery, New York, NY, USA, 1–12. https://doi.org/10.1145/3173574.3173735[13] Thomas Langerak, Juan José Zárate, David Lindlbauer, Christian Holz, and Otmar
Hilliges. 2020. Omni: Volumetric Sensing and Actuation of Passive Magnetic
Tools for Dynamic Haptic Feedback. In Proceedings of the 33rd Annual ACM
Symposium on User Interface Software and Technology (Virtual Event, USA) (UIST’20). Association for Computing Machinery, New York, NY, USA, 594–606. https://doi.org/10.1145/3379337.3415589[14] Karon E MacLean. 2008. Haptic interaction design for everyday interfaces.
Reviews of Human Factors and Ergonomics 4, 1 (2008), 149–194.
[15] Luigi Malomo, Jesús Pérez, Emmanuel Iarussi, Nico Pietroni, Eder Miguel, Paolo
Cignoni, and Bernd Bickel. 2018. FlexMaps: Computational Design of Flat Flexible
Shells for Shaping 3D Objects. ACM Trans. Graph. 37, 6, Article 241 (Dec. 2018),
14 pages. https://doi.org/10.1145/3272127.3275076[16] Xingyu Ni, Bo Zhu, Bin Wang, and Baoquan Chen. 2020. A Level-Set Method forMagnetic Substance Simulation. ACM Trans. Graph. 39, 4, Article 29 (July 2020),
15 pages. https://doi.org/10.1145/3386569.3392445[17] Jorge Nocedal and Stephen J. Wright. 2006. Numerical Optimization (2nd ed.).
Springer Science+Business Media. https://doi.org/10.1007/978-0-387-40065-5[18] Masa Ogata. 2018. Magneto-Haptics: Embedding Magnetic Force Feedback forPhysical Interactions. In Proceedings of the 31st Annual ACM Symposium on User
Interface Software and Technology (Berlin, Germany) (UIST ’18). Association for
Computing Machinery, New York, NY, USA, 737–743. https://doi.org/10.1145/
3242587.3242615[19] Masa Ogata and Masaaki Fukumoto. 2015. FluxPaper: Reinventing Paper with
Dynamic Actuation Powered by Magnetic Flux. In Proceedings of the 33rd Annual
ACM Conference on Human Factors in Computing Systems (Seoul, Republic ofKorea) (CHI ’15). ACM, New York, NY, USA, 29–38. https://doi.org/10.1145/
2702123.2702516[20] Antti Oulasvirta, Niraj Ramesh Dayama, Morteza Shiripour, Maximilian John,
and Andreas Karrenbauer. 2020. Combinatorial Optimization of Graphical User
Interface Designs. Proc. IEEE 108, 3 (March 2020), 434–464. https://doi.org/10.
1109/JPROC.2020.2969687[21] Jia Pan, Sachin Chitta, and Dinesh Manocha. 2012. FCL: A General Purpose Library for Collision and Proximity Queries. In 2012 IEEE International Conference on
Robotics and Automation. 3859–3866. https://doi.org/10.1109/ICRA.2012.6225337[22] Michal Piovarči, Danny M. Kaufman, David I. W. Levin, and Piotr Didyk. 2020.
Fabrication-in-the-Loop Co-Optimization of Surfaces and Styli for Drawing
Haptics. ACM Trans. Graph. 39, 4, Article 116 (July 2020), 16 pages. https://doi.org/10.1145/3386569.3392467[23] Michal Piovarči, David I. W. Levin, Danny M. Kaufman, and Piotr Didyk. 2018.
Perception-Aware Modeling and Fabrication of Digital Drawing Tools. ACM

CHI ’21, May 8–13, 2021, Yokohama, Japan

Trans. Graph. 37, 4, Article 123 (July 2018), 15 pages. https://doi.org/10.1145/
3197517.3201322[24] Michal Piovarči, David I. W. Levin, Jason Rebello, Desai Chen, Roman Ďurikovič,
Hanspeter Pfister, Wojciech Matusik, and Piotr Didyk. 2016. An InteractionAware, Perceptual Model for Non-Linear Elastic Objects. ACM Trans. Graph. 35,
4, Article 55 (July 2016), 13 pages. https://doi.org/10.1145/2897824.2925885[25] Romain Prévost, Emily Whiting, Sylvain Lefebvre, and Olga Sorkine-Hornung.
2013. Make It Stand: Balancing Shapes for 3D Fabrication. ACM Trans. Graph.
32, 4, Article 81 (July 2013), 10 pages. https://doi.org/10.1145/2461912.2461957[26] Gabriel Robles-De-La-Torre and Vincent Hayward. 2001. Force can overcome
object geometry in the perception of shape through active touch. Nature 412,
6845 (01 Jul 2001), 445–448. https://doi.org/10.1038/35086588[27] Christian Schumacher, Bernd Bickel, Jan Rys, Steve Marschner, Chiara Daraio,
and Markus Gross. 2015. Microstructures to Control Elasticity in 3D Printing.
ACM Trans. Graph. 34, 4, Article 136 (July 2015), 13 pages. https://doi.org/10.
1145/2766926[28] Bobak Shahriari, Kevin Swersky, Ziyu Wang, Ryan P. Adams, and Nando deFreitas. 2016. Taking the Human Out of the Loop: A Review of Bayesian Optimization. Proc. IEEE 104, 1 (January 2016), 148–175. https://doi.org/10.1109/
JPROC.2015.2494218[29] Mélina Skouras, Bernhard Thomaszewski, Stelian Coros, Bernd Bickel, andMarkus Gross. 2013. Computational Design of Actuated Deformable Characters. ACM Trans. Graph. 32, 4, Article 82 (July 2013), 10 pages. https://doi.org/10.1145/2461912.2461979[30] Bernhard Thomaszewski, Andreas Gumann, Simon Pabst, and Wolfgang Straßer.
2008. Magnets in Motion. ACM Trans. Graph. 27, 5, Article 162 (Dec. 2008),
9 pages. https://doi.org/10.1145/1409060.1409115[31] Chelsea Tymms, Esther P. Gardner, and Denis Zorin. 2018. A Quantitative
Perceptual Model for Tactile Roughness. ACM Trans. Graph. 37, 5, Article 168(Oct. 2018), 14 pages. https://doi.org/10.1145/3186267[32] Nobuyuki Umetani, Yuki Koyama, Ryan Schmidt, and Takeo Igarashi. 2014.
Pteromys: Interactive Design and Optimization of Free-Formed Free-Flight
Model Airplanes. ACM Trans. Graph. 33, 4, Article 65 (July 2014), 10 pages.
https://doi.org/10.1145/2601097.2601129[33] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler
Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser,
Jonathan Bright, Stéfan J. van der Walt, Matthew Brett, Joshua Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Eric Jones, Robert Kern,
Eric Larson, CJ Carey, İlhan Polat, Yu Feng, Eric W. Moore, Jake Vand erPlas,
Denis Laxalde, Josef Perktold, Robert Cimrman, Ian Henriksen, E. A. Quintero,Charles R Harris, Anne M. Archibald, Antônio H. Ribeiro, Fabian Pedregosa, Paulvan Mulbregt, and SciPy 1.0 Contributors. 2020. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods 17 (2020), 261–272.
https://doi.org/10.1038/s41592-019-0686-2[34] Malte Weiss, Chat Wacharamanotham, Simon Voelker, and Jan Borchers. 2011.
FingerFlux: Near-surface Haptic Feedback on Tabletops. In Proceedings of the24th Annual ACM Symposium on User Interface Software and Technology (Santa
Barbara, California, USA) (UIST ’11). ACM, New York, NY, USA, 615–620. https://doi.org/10.1145/2047196.2047277[35] Kentaro Yasu. 2017. Magnetic Plotter: A Macrotexture Design Method Using
Magnetic Rubber Sheets. In Proceedings of the 2017 CHI Conference on Human
Factors in Computing Systems (Denver, Colorado, USA) (CHI ’17). Association for
Computing Machinery, New York, NY, USA, 4983–4993. https://doi.org/10.1145/
3025453.3025702[36] Kentaro Yasu. 2020. MagneLayer: Force Field Fabrication by Layered Magnetic
Sheets. In Proceedings of the 2020 CHI Conference on Human Factors in Computing
Systems (Honolulu, HI, USA) (CHI ’20). Association for Computing Machinery,
New York, NY, USA, 1–9. https://doi.org/10.1145/3313831.3376552[37] Kar W Yung, Peter B Landecker, and Daniel D Villani. 1998. An analytic solution
for the force between two magnetic dipoles. Physical Separation in Science andEngineering 9, 1 (1998), 39–52.
[38] Jonas Zehnder, Espen Knoop, Moritz Bächer, and Bernhard Thomaszewski. 2017.
Metasilicone: Design and Fabrication of Composite Silicone with Desired Mechanical Properties. ACM Trans. Graph. 36, 6, Article 240 (Nov. 2017), 13 pages.
https://doi.org/10.1145/3130800.3130881[39] Clement Zheng and Ellen Yi-Luen Do. 2018. Mechamagnets: Tactile Mechanisms
with Embedded Magnets. In Proceedings of the Twelfth International Conferenceon Tangible, Embedded, and Embodied Interaction (Stockholm, Sweden) (TEI ’18).
ACM, New York, NY, USA, 57–64. https://doi.org/10.1145/3173225.3173268