mechanicalForcesUdi.cpp
59f89f10a1ebmaster

mechanicalForcesUdi.cpp
View Options

			#include "mechanicalForcesUdi.hpp"

			using namespace std;


			void MechanicalForcesUdi::operator()(const state_type &y,state_type &dydt,const double t){
			fromState(y);
			fill(dydt.begin(), dydt.end(), 0);
			uint nbcells = lattice->getNbCells();
			for (uint c=0;c<nbcells;c++){
			Cell& ce = lattice->getCell(c);
			dSpecialEnergy(&ce,y,dydt);
			}
			// double n2= norm2(dydt);
			// not sure if needed...
			//checkConvexity(y,dydt);

			}

			void MechanicalForcesUdi::dSpecialEnergy(Cell *c,const state_type &y,state_type &dydt){
			BoundaryIterator it(c);
			Vec2 from1,to1;
			float darea = area(c) - area0(c);
			for(Boundary *b=it.first();!it.isLast();b=it.next()){
			uint fidx = b->from()->vertexStateIndex();
			uint tidx = b->to()->vertexStateIndex();
			b->relativePosition(y[fidx],y[fidx+1],y[tidx],y[tidx+1],from1, to1);
			Vec2 v= to1-from1;
			if(v.Norm()>0.0001){
			v.Normalize();
			dydt[tidx] -= (k_perim * c->perimeter() + 0.5k_areadarea) * v[0];
			dydt[tidx+1] -= (k_perim * c->perimeter() + 0.5k_areadarea) * v[0];
			}
			else{
			dydt[tidx] -= 1;
			dydt[tidx+1] -= 1;
			}
			}
			}

			void MechanicalForcesUdi::updateSignal(const state_type &y){}