ComputeSphericalBC.C

Go to the documentation of this file.

#include "ComputeSphericalBC.h"
#include "Node.h"
#include "SimParameters.h"
#include "Patch.h"

/************************************************************************/
/*                                                                      */
/*                      FUNCTION ComputeSphericalBC                     */
/*                                                                      */
/*      This is the constructor for the ComputeSphericalBC force object.*/
/*   It is responsible for getting all the parameters from the          */
/*   SimParameters object and then determining if this object needs     */
/*   to perform any computation.  It only needs to do so if there is    */
/*   some portion of the patch that lays outside of the spherical       */
/*   boundaries.                                                        */
/*                                                                      */
/************************************************************************/

ComputeSphericalBC::ComputeSphericalBC(ComputeID c, PatchID pid)
  : ComputeHomePatch(c,pid)
{
        reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);

        SimParameters *simParams = Node::Object()->simParameters;

        //  Get parameters from the SimParameters object
        r1 = simParams->sphericalBCr1;
        r2 = simParams->sphericalBCr2;
        r1_2 = r1*r1;
        r2_2 = r2*r2;
        k1 = simParams->sphericalBCk1;
        k2 = simParams->sphericalBCk2;
        exp1 = simParams->sphericalBCexp1;
        exp2 = simParams->sphericalBCexp2;

        //  Check to see if this is one set of parameters or two
        if (r2 > -1.0)
        {
                twoForces = TRUE;
        }
        else
        {
                twoForces = FALSE;
        }

        center = simParams->sphericalCenter;

}
/*                      END OF FUNCTION ComputeSphericalBC              */

/************************************************************************/
/*                                                                      */
/*                      FUNCTION ~ComputeSphericalBC                    */
/*                                                                      */
/*      This is the destructor for the ComputeSphericalBC force object. */
/*   It currently does ABSOLUTELY NOTHING!!                             */
/*                                                                      */
/************************************************************************/

ComputeSphericalBC::~ComputeSphericalBC()

{
        delete reduction;
}
/*                      END OF FUNCTION ~ComputeSphericalBC             */

/************************************************************************/
/*                                                                      */
/*                              FUNCTION force                          */
/*                                                                      */
/*   INPUTS:                                                            */
/*      numAtoms - Number of coordinates being passed                   */
/*      x - Array of atom coordinates                                   */
/*      forces - Array of force vectors                                 */
/*                                                                      */
/*      This function calculates the force and energy for the shereical */
/*   boundary conditions for this patch.                                */
/*                                                                      */
/************************************************************************/
void ComputeSphericalBC::doForce(FullAtom* p, Results* r)
{
        Vector diff;            //  Distance from atom to center of sphere
        Vector f;               //  Calculated force vector
        int i, j;               //  Loop counters
        BigReal dist, dist_2;   //  Distance from atom to center, and this
                                //  distance squared
        BigReal rval;           //  Difference between distance from atom
                                //  to center and radius of sphere
        BigReal eval;           //  Energy value for this atom
        BigReal fval;           //  Force magnitude for this atom

        // aliases to work with old code
        FullAtom *x = p;
        Force *forces = r->f[Results::normal];
        BigReal energy = 0;

        //  Loop through and check each atom
        for (i=0; i<numAtoms; i++)
        {
                //  Calculate the vector from the atom to the center of the
                //  sphere
                diff.x = x[i].position.x - center.x;
                diff.y = x[i].position.y - center.y;
                diff.z = x[i].position.z - center.z;
                
                //  Calculate the distance squared
                dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;

                //  Look to see if we are outside either radius
                if ( (dist_2 > r1_2) || (twoForces && (dist_2 > r2_2)) )
                {
                        //  Calculate the distance to the center
                        dist = sqrt(dist_2);

                        //  Normalize the direction vector
                        diff /= dist;

                        //  Check to see if we are outside radius 1
                        if (dist > r1)
                        {
                                //  Assign the force vector to the
                                //  unit direction vector
                                f.x = diff.x;
                                f.y = diff.y;
                                f.z = diff.z;

                                //  Calculate the energy which is
                                //  e = k1*(r_i-r_center)^exp1
                                eval = k1;
                                rval = fabs(dist - r1);

                                for (j=0; j<exp1; j++)
                                {
                                        eval *= rval;
                                }

                                energy += eval;

                                //  Now calculate the force which is
                                //  e = -k1*exp1*(r_i-r_center1)^(exp1-1)
                                fval = -exp1*k1;

                                for (j=0; j<exp1-1; j++)
                                {
                                        fval *= rval;
                                }

                                //  Multiply the force magnitude to the
                                //  unit direction vector to get the
                                //  resulting force
                                f *= fval;

                                //  Add the force to the force vectors
                                forces[i].x += f.x;
                                forces[i].y += f.y;
                                forces[i].z += f.z;
                        }

                        //  Check to see if there is a second radius
                        //  and if we are outside of it
                        if (twoForces && (dist > r2) )
                        {
                                //  Assign the force vector to the
                                //  unit direction vector
                                f.x = diff.x;
                                f.y = diff.y;
                                f.z = diff.z;

                                //  Calculate the energy which is
                                //  e = k2*(r_i-r_center2)^exp2
                                eval = k2;
                                rval = fabs(dist - r2);

                                for (j=0; j<exp2; j++)
                                {
                                        eval *= rval;
                                }

                                energy += eval;

                                //  Now calculate the force which is
                                //  e = -k2*exp2*(r_i-r_center2)^(exp2-1)
                                fval = -exp2*k2;

                                for (j=0; j<exp2-1; j++)
                                {
                                        fval *= rval;
                                }

                                //  Multiply the force magnitude to the
                                //  unit direction vector to get the
                                //  resulting force
                                f *= fval;

                                //  Add the force to the force vectors
                                forces[i].x += f.x;
                                forces[i].y += f.y;
                                forces[i].z += f.z;
                        }
                }
        }

    reduction->item(REDUCTION_BC_ENERGY) += energy;
    reduction->submit();

}
/*                      END OF FUNCTION force                           */