Commit a89be473 by Francois Gygi

Added counting of blocked dofs in constraints


git-svn-id: http://qboxcode.org/svn/qb/trunk@698 cba15fb0-1239-40c8-b417-11db7ca47a34
parent ff6b74e2
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// AngleConstraint.h // AngleConstraint.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: AngleConstraint.h,v 1.7 2008-11-14 04:01:25 fgygi Exp $ // $Id: AngleConstraint.h,v 1.8 2009-05-15 04:38:48 fgygi Exp $
#ifndef ANGLECONSTRAINT_H #ifndef ANGLECONSTRAINT_H
#define ANGLECONSTRAINT_H #define ANGLECONSTRAINT_H
...@@ -75,6 +75,7 @@ class AngleConstraint : public Constraint ...@@ -75,6 +75,7 @@ class AngleConstraint : public Constraint
} }
void setup(const AtomSet& atoms); void setup(const AtomSet& atoms);
int dofs(void) const { return 1; }
void update(double dt); void update(double dt);
bool enforce_r(const std::vector<std::vector<double> > &r0, bool enforce_r(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &rp) const; std::vector<std::vector<double> > &rp) const;
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// Constraint.h // Constraint.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: Constraint.h,v 1.7 2008-11-14 04:01:26 fgygi Exp $ // $Id: Constraint.h,v 1.8 2009-05-15 04:38:48 fgygi Exp $
#ifndef CONSTRAINT_H #ifndef CONSTRAINT_H
#define CONSTRAINT_H #define CONSTRAINT_H
...@@ -52,6 +52,7 @@ class Constraint ...@@ -52,6 +52,7 @@ class Constraint
const std::vector<std::vector<double> > &f) = 0; const std::vector<std::vector<double> > &f) = 0;
virtual void update(double dt) = 0; virtual void update(double dt) = 0;
virtual void setup(const AtomSet& atoms) = 0; virtual void setup(const AtomSet& atoms) = 0;
virtual int dofs(void) const = 0;
virtual std::ostream& print(std::ostream &os) = 0; virtual std::ostream& print(std::ostream &os) = 0;
std::string name(void) const { return name_; } std::string name(void) const { return name_; }
std::string names(int i) const std::string names(int i) const
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// ConstraintSet.C // ConstraintSet.C
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintSet.C,v 1.8 2009-04-30 22:24:51 fgygi Exp $ // $Id: ConstraintSet.C,v 1.9 2009-05-15 04:38:48 fgygi Exp $
#include "ConstraintSet.h" #include "ConstraintSet.h"
#include "PositionConstraint.h" #include "PositionConstraint.h"
...@@ -435,6 +435,9 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv) ...@@ -435,6 +435,9 @@ bool ConstraintSet::define_constraint(AtomSet &atoms, int argc, char **argv)
return false; return false;
} }
// update total number of blocked degrees of freedom
ndofs_ += constraint_list.back()->dofs();
return true; return true;
} }
...@@ -505,6 +508,10 @@ bool ConstraintSet::delete_constraint(int argc, char **argv) ...@@ -505,6 +508,10 @@ bool ConstraintSet::delete_constraint(int argc, char **argv)
if ( pc->name() == name ) if ( pc->name() == name )
{ {
found = true; found = true;
// update total number of blocked degrees of freedom
ndofs_ -= pc->dofs();
delete pc; delete pc;
// remove constraint pointer from the list // remove constraint pointer from the list
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// ConstraintSet.h // ConstraintSet.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: ConstraintSet.h,v 1.7 2008-11-14 04:01:26 fgygi Exp $ // $Id: ConstraintSet.h,v 1.8 2009-05-15 04:38:48 fgygi Exp $
#ifndef CONSTRAINTSET_H #ifndef CONSTRAINTSET_H
#define CONSTRAINTSET_H #define CONSTRAINTSET_H
...@@ -34,6 +34,8 @@ class ConstraintSet ...@@ -34,6 +34,8 @@ class ConstraintSet
const Context& ctxt_; const Context& ctxt_;
std::vector<Constraint *> constraint_list; std::vector<Constraint *> constraint_list;
// ndofs_: total number of degrees of freedom blocked by the constraints
int ndofs_;
public: public:
...@@ -44,6 +46,7 @@ class ConstraintSet ...@@ -44,6 +46,7 @@ class ConstraintSet
bool delete_constraint(int argc, char **argv); bool delete_constraint(int argc, char **argv);
void list_constraints(std::ostream &os); void list_constraints(std::ostream &os);
int size(void) const { return constraint_list.size(); } int size(void) const { return constraint_list.size(); }
int ndofs(void) const { return ndofs_; }
void enforce(AtomSet& atoms); void enforce(AtomSet& atoms);
void enforce_r(const std::vector<std::vector<double> > &r0, void enforce_r(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &rp); std::vector<std::vector<double> > &rp);
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// DistanceConstraint.h // DistanceConstraint.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: DistanceConstraint.h,v 1.7 2008-11-14 04:01:26 fgygi Exp $ // $Id: DistanceConstraint.h,v 1.8 2009-05-15 04:38:48 fgygi Exp $
#ifndef DISTANCECONSTRAINT_H #ifndef DISTANCECONSTRAINT_H
#define DISTANCECONSTRAINT_H #define DISTANCECONSTRAINT_H
...@@ -65,6 +65,7 @@ class DistanceConstraint : public Constraint ...@@ -65,6 +65,7 @@ class DistanceConstraint : public Constraint
} }
void setup(const AtomSet& atoms); void setup(const AtomSet& atoms);
int dofs(void) const { return 1; }
void update(double dt); void update(double dt);
bool enforce_r(const std::vector<std::vector<double> > &r0, bool enforce_r(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &rp) const; std::vector<std::vector<double> > &rp) const;
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// IonicStepper.h: // IonicStepper.h:
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: IonicStepper.h,v 1.11 2008-09-08 15:56:18 fgygi Exp $ // $Id: IonicStepper.h,v 1.12 2009-05-15 04:39:07 fgygi Exp $
#ifndef IONICSTEPPER_H #ifndef IONICSTEPPER_H
#define IONICSTEPPER_H #define IONICSTEPPER_H
...@@ -33,7 +33,8 @@ class IonicStepper ...@@ -33,7 +33,8 @@ class IonicStepper
ConstraintSet& constraints_; ConstraintSet& constraints_;
double dt_; double dt_;
int nsp_; int nsp_;
int ndofs_; // ndofs_ is the total number of degrees of freedom after constraints are considered
int ndofs_;
std::vector<int> na_; // number of atoms per species na_[nsp_] std::vector<int> na_; // number of atoms per species na_[nsp_]
std::vector<std::vector< double> > r0_; // r0_[nsp_][3*na_] std::vector<std::vector< double> > r0_; // r0_[nsp_][3*na_]
std::vector<std::vector< double> > rp_; // rp_[nsp_][3*na_] std::vector<std::vector< double> > rp_; // rp_[nsp_][3*na_]
...@@ -46,7 +47,7 @@ class IonicStepper ...@@ -46,7 +47,7 @@ class IonicStepper
IonicStepper (Sample& s) : s_(s), atoms_(s.atoms), IonicStepper (Sample& s) : s_(s), atoms_(s.atoms),
constraints_(s.constraints), dt_(s.ctrl.dt) constraints_(s.constraints), dt_(s.ctrl.dt)
{ {
ndofs_ = 3 * atoms_.size() - constraints_.size(); ndofs_ = 3 * atoms_.size() - constraints_.ndofs();
// if there are more constraints than dofs, set ndofs_ to zero // if there are more constraints than dofs, set ndofs_ to zero
if ( ndofs_ < 0 ) ndofs_ = 0; if ( ndofs_ < 0 ) ndofs_ = 0;
nsp_ = atoms_.nsp(); nsp_ = atoms_.nsp();
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// PositionConstraint.h // PositionConstraint.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: PositionConstraint.h,v 1.1 2009-04-30 23:13:39 fgygi Exp $ // $Id: PositionConstraint.h,v 1.2 2009-05-15 04:38:48 fgygi Exp $
#ifndef POSITIONCONSTRAINT_H #ifndef POSITIONCONSTRAINT_H
#define POSITIONCONSTRAINT_H #define POSITIONCONSTRAINT_H
...@@ -61,6 +61,7 @@ class PositionConstraint : public Constraint ...@@ -61,6 +61,7 @@ class PositionConstraint : public Constraint
} }
void setup(const AtomSet& atoms); void setup(const AtomSet& atoms);
int dofs(void) const { return 3; }
void update(double dt); void update(double dt);
bool enforce_r(const std::vector<std::vector<double> > &r0, bool enforce_r(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &rp) const; std::vector<std::vector<double> > &rp) const;
......
...@@ -15,7 +15,7 @@ ...@@ -15,7 +15,7 @@
// TorsionConstraint.h // TorsionConstraint.h
// //
//////////////////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////
// $Id: TorsionConstraint.h,v 1.7 2008-11-14 04:01:26 fgygi Exp $ // $Id: TorsionConstraint.h,v 1.8 2009-05-15 04:38:48 fgygi Exp $
#ifndef TORSIONCONSTRAINT_H #ifndef TORSIONCONSTRAINT_H
#define TORSIONCONSTRAINT_H #define TORSIONCONSTRAINT_H
...@@ -80,6 +80,7 @@ class TorsionConstraint : public Constraint ...@@ -80,6 +80,7 @@ class TorsionConstraint : public Constraint
} }
void setup(const AtomSet& atoms); void setup(const AtomSet& atoms);
int dofs(void) const { return 1; }
void update(double dt); void update(double dt);
bool enforce_r(const std::vector<std::vector<double> > &r0, bool enforce_r(const std::vector<std::vector<double> > &r0,
std::vector<std::vector<double> > &rp) const; std::vector<std::vector<double> > &rp) const;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment