The force class have a pre-determined interface, set by an abstract base class, to allow them to always be called from fix_cpl_force in the same way. Think of an interface like a contract, you guarentee that you will always take the same arguments in and return the same things. In this case, we have an interface which takes in position, velocity, acceleration, mass, radius (sigma) and interaction strength (epsilon) and works out summed up arrays needed for the particular force type (in the pre_force function) and returns the actual force (in the get_force function). As each type of coupling force you could want to apply always has the same form:

0) Constructor : Create the force class with all required argments stored in a std::map. 1) pre_force : Get some stuff from the particle system (e.g. add up current void fraction, get latest velocity from CFD solver) 2) get_force : Use the pre_force information, with the CFD solver information, to get the force on each particle.

These means we can make use of c++ polymorphism, where we choose which type of force we want based on the forcetype input argument. The required force type object is instantiated using a "factory" which takes the user input and returns fxyz, a pointer to the appropriate object: 

//Force factory
if (fxyzType.compare("Flekkoy") == 0) {
    fxyz.reset(new CPLForceFlekkoy(9, cfdBuf->shape(1), 
                                      cfdBuf->shape(2), 
                                      cfdBuf->shape(3)));
} else if (fxyzType.compare("test") == 0) {
    fxyz.reset(new CPLForceTest(3, cfdBuf->shape(1), 
                                   cfdBuf->shape(2), 
                                   cfdBuf->shape(3)));
} else if (fxyzType.compare("Velocity") == 0) {
    fxyz.reset(new CPLForceVelocity(3, cfdBuf->shape(1), 
                                       cfdBuf->shape(2), 
                                       cfdBuf->shape(3)));
} else if (fxyzType.compare("Drag") == 0) {
    fxyz.reset(new CPLForceDrag(9, cfdBuf->shape(1), 
                                   cfdBuf->shape(2), 
                                   cfdBuf->shape(3),
                                   arg_map)); 

} else if (fxyzType.compare("Di_Felice") == 0) {
    fxyz.reset(new CPLForceGranular(9, cfdBuf->shape(1), 
                                       cfdBuf->shape(2), 
                                       cfdBuf->shape(3),
                                       arg_map)); 

} else {
    std::string cmd("CPLForce type ");
    cmd += fxyzType + " not defined";
    throw std::runtime_error(cmd);
}

The arg_map is a std::map, basically a set of paired {keywords and values}. These are obtained directly from the user input taking alternating keywords and values following forcetype. For example, if you had specified an input of the form: 

fix 5 all cpl/init region all forcetype Di_Felice Cd 0.0005 overlap false me 1e-4 rho 1e3 sendtype velocity

then arg_map would be built up by parsing the commands after "Di_Felice" as 

 Key           value
"Cd"            "0.0005"
"overlap"        "false"
"mu"            "1e-4"
"rho"           "1e3"

and this would be passed to the Di_Felice force type, instantiating a CPLForceGranular force object and setting the pointer fxyz to this. Once we have this fxyz pointer, it can then be used by looping over all particles in LAMMPS (nlocal here) 

    //Pre-force calculation, get quantities from discrete system needed to apply force
	for (int i = 0; i < nlocal; ++i)
	{
   		if (mask[i] & groupbit)
    	{
	        //Get local molecule data
	        mi = rmass[i];
	        radi = radius[i];
	        for (int n=0; n<3; n++){
	            xi[n]=x[i][n]; 
	            vi[n]=v[i][n]; 
	            ai[n]=f[i][n];
	        }

	        // Sum all the weights for each cell.
	        fxyz->pre_force(xi, vi, ai, mi, radi, pot);

    	}
    }


//This pre-force has added all the needed  things for the force you want, so we can simply get the force now:

   // Calculate force and apply
    for (int i = 0; i < nlocal; ++i)
    {
        if (mask[i] & groupbit)
        {

            //Get local molecule data
            mi = rmass[i];
            radi = radius[i];
            for (int n=0; n<3; n++){
                xi[n]=x[i][n]; 
                vi[n]=v[i][n]; 
                ai[n]=f[i][n];
            }

            //Get force from object
            fi = fxyz->get_force(xi, vi, ai, mi, radi, pot);

        }
    }