[1] This is the first paper in a two-part series examining the fluid dynamics of crystal settling and flotation in the lunar magma ocean. A key challenge in constraining solidification processes is determining the ability of individual crystals to decouple from vigorous thermal convection and settle out or float. The goal of this paper is to develop a computational methodology capable of capturing the complex solid-fluid interactions that determine settling and flotation. In the second paper, we use this computational approach to explore the conditions under which plagioclase feldspar would be able to buoyantly float and form the earliest crust on the Moon. The direct numerical method described in this paper relies on a fictitious domain approach and captures solid-body motion in 2D and 3D with little overhead beyond single fluid calculations. The two main innovations of our numerical implementation of a fictitious domain approach are an analytical quadrature scheme, which increases accuracy and reduces computational expense, and the derivation of a multibody collision scheme. Advantages of this approach over previous simulations of crystal-bearing magmatic suspensions include the following: (1) we fully resolve the two-way interaction between fluid and solid phases, implying that crystals are not only passively advected in an ambient flow field but are also actively driving flow, and (2) we resolve the flow around each individual crystal without assuming specific settling speeds or drag coefficients. We present several benchmark problems and convergence tests to validate our approach.