src/integrators/NVE.cpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */
 
 /**
  * @file NVE.cpp
  * @author tlin
  * @date 11/08/2004
  * @time 15:13am
  * @version 1.0
  */

#include "integrators/IntegratorCreator.hpp"  
#include "integrators/NVE.hpp"
#include "primitives/Molecule.hpp"
#include "utils/OOPSEConstant.hpp"

namespace oopse {

static IntegratorBuilder<NVE>* NVECreator = new IntegratorBuilder<NVE>("NVE");

NVE::NVE(SimInfo* info) : VelocityVerletIntegrator(info){

}

void NVE::moveA(){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    Vector3d vel;
    Vector3d pos;
    Vector3d frc;
    Vector3d Tb;
    Vector3d ji;
    double mass;
    
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {

                vel =integrableObject->getVel();
                pos = integrableObject->getPos();
                frc = integrableObject->getFrc();
                mass = integrableObject->getMass();
                
                // velocity half step
                vel += (dt2 /mass * OOPSEConstant::energyConvert) * frc;

                // position whole step
                pos += dt * vel;

                integrableObject->setVel(vel);
                integrableObject->setPos(pos);

                if (integrableObject->isDirectional()){

                    // get and convert the torque to body frame

                    Tb = integrableObject->lab2Body(integrableObject->getTrq());

                    // get the angular momentum, and propagate a half step

                    ji = integrableObject->getJ();

                    ji += (dt2  * OOPSEConstant::energyConvert) * Tb;

                    rotAlgo->rotate(integrableObject, ji, dt);

                    integrableObject->setJ(ji);
                }

            
        }
    } //end for(mol = info_->beginMolecule(i))
    
    //constraintAlgorithm->doConstrainA();
    
}    

void NVE::moveB(){
    SimInfo::MoleculeIterator i;
    Molecule::IntegrableObjectIterator  j;
    Molecule* mol;
    StuntDouble* integrableObject;
    Vector3d vel;
    Vector3d frc;
    Vector3d Tb;
    Vector3d ji;
    double mass;
    
    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
               integrableObject = mol->nextIntegrableObject(j)) {

                vel =integrableObject->getVel();
                frc = integrableObject->getFrc();
                mass = integrableObject->getMass();
                
                // velocity half step
                vel += (dt2 /mass * OOPSEConstant::energyConvert) * frc;
                
                integrableObject->setVel(vel);

                if (integrableObject->isDirectional()){

                    // get and convert the torque to body frame

                    Tb = integrableObject->lab2Body(integrableObject->getTrq());

                    // get the angular momentum, and propagate a half step

                    ji = integrableObject->getJ();

                    ji += (dt2  * OOPSEConstant::energyConvert) * Tb;

                    integrableObject->setJ(ji);
                }

            
        }
    } //end for(mol = info_->beginMolecule(i))
  

    //constraintAlgorithm->doConstrainB();

}


double NVE::calcConservedQuantity() {
    return thermo.getTotalE();
}

} //end namespace oopse
Revision:	1837
Committed:	Thu Dec 2 22:15:31 2004 UTC (19 years, 7 months ago) by tim
File size:	4736 byte(s)
Log Message:	refine factory pattern to make it initialized correctly
#	User	Rev	Content
1	tim	1722	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			/**
27			* @file NVE.cpp
28			* @author tlin
29			* @date 11/08/2004
30			* @time 15:13am
31			* @version 1.0
32			*/
33	tim	1837
34			#include "integrators/IntegratorCreator.hpp"
35	tim	1820	#include "integrators/NVE.hpp"
36			#include "primitives/Molecule.hpp"
37			#include "utils/OOPSEConstant.hpp"
38	tim	1722
39			namespace oopse {
40
41	tim	1837	static IntegratorBuilder<NVE>* NVECreator = new IntegratorBuilder<NVE>("NVE");
42	tim	1758
43	tim	1756	NVE::NVE(SimInfo* info) : VelocityVerletIntegrator(info){
44	tim	1722
45			}
46
47	tim	1756	void NVE::moveA(){
48	tim	1820	SimInfo::MoleculeIterator i;
49			Molecule::IntegrableObjectIterator j;
50	tim	1722	Molecule* mol;
51			StuntDouble* integrableObject;
52			Vector3d vel;
53			Vector3d pos;
54			Vector3d frc;
55			Vector3d Tb;
56			Vector3d ji;
57			double mass;
58
59	tim	1756	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
60	tim	1722	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
61			integrableObject = mol->nextIntegrableObject(j)) {
62
63			vel =integrableObject->getVel();
64			pos = integrableObject->getPos();
65			frc = integrableObject->getFrc();
66			mass = integrableObject->getMass();
67
68			// velocity half step
69	tim	1820	vel += (dt2 /mass * OOPSEConstant::energyConvert) * frc;
70	tim	1722
71			// position whole step
72	tim	1774	pos += dt * vel;
73	tim	1722
74			integrableObject->setVel(vel);
75			integrableObject->setPos(pos);
76
77			if (integrableObject->isDirectional()){
78
79			// get and convert the torque to body frame
80
81			Tb = integrableObject->lab2Body(integrableObject->getTrq());
82
83			// get the angular momentum, and propagate a half step
84
85			ji = integrableObject->getJ();
86
87	tim	1820	ji += (dt2 * OOPSEConstant::energyConvert) * Tb;
88	tim	1722
89	tim	1820	rotAlgo->rotate(integrableObject, ji, dt);
90	tim	1722
91			integrableObject->setJ(ji);
92			}
93
94
95			}
96	tim	1756	} //end for(mol = info_->beginMolecule(i))
97	tim	1722
98	tim	1758	//constraintAlgorithm->doConstrainA();
99	tim	1722
100			}
101
102	tim	1756	void NVE::moveB(){
103	tim	1820	SimInfo::MoleculeIterator i;
104			Molecule::IntegrableObjectIterator j;
105	tim	1722	Molecule* mol;
106			StuntDouble* integrableObject;
107			Vector3d vel;
108			Vector3d frc;
109			Vector3d Tb;
110			Vector3d ji;
111			double mass;
112
113	tim	1756	for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
114	tim	1722	for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
115			integrableObject = mol->nextIntegrableObject(j)) {
116
117			vel =integrableObject->getVel();
118			frc = integrableObject->getFrc();
119			mass = integrableObject->getMass();
120
121			// velocity half step
122	tim	1820	vel += (dt2 /mass * OOPSEConstant::energyConvert) * frc;
123	tim	1722
124			integrableObject->setVel(vel);
125
126			if (integrableObject->isDirectional()){
127
128			// get and convert the torque to body frame
129
130			Tb = integrableObject->lab2Body(integrableObject->getTrq());
131
132			// get the angular momentum, and propagate a half step
133
134			ji = integrableObject->getJ();
135
136	tim	1820	ji += (dt2 * OOPSEConstant::energyConvert) * Tb;
137	tim	1722
138			integrableObject->setJ(ji);
139			}
140
141
142			}
143	tim	1756	} //end for(mol = info_->beginMolecule(i))
144	tim	1722
145
146	tim	1758	//constraintAlgorithm->doConstrainB();
147	tim	1722
148			}
149
150	tim	1756
151	tim	1774	double NVE::calcConservedQuantity() {
152	tim	1820	return thermo.getTotalE();
153	tim	1756	}
154
155	tim	1722	} //end namespace oopse