In this paper, we present a parallel computational tool, BeamBeam3D, developed at Lawrence Berkeley National Laboratory, for strong-strong/strong-weak beam-beam modeling. This tool calculates self-consistently the electromagnetic beam-beam forces for arbitrary distributions during each collision when a strong-strong beam-beam interaction model is used. When a strong-weak model is used, the code has the option of using a Gaussian approximation for the strong beam. BeamBeam3D uses a multiple-slice model, so finite bunch length effects can be studied. The code also includes a Lorentz boost and rotation to treat collisions with finite collision crossing angle. It handles arbitrary closed-orbit separation (static or time dependent) and models long-range beam-beam interactions using a newly developed shifted Green function approach. It can also handle multiple interaction points using externally supplied linear maps between interaction points in the strong-weak model. The code has been used to study beam-beam effects in the RHIC, Tevatron, and LHC. In this paper we will describe the BeamBeam3D code, present example simulations, and describe the code performance.