moving.cpp_bak

#include "ISOP2P1.h"
#include "preconditioner.h"
#include "functions.h"

#define DIM 2

void ISOP2P1::syncMesh()
{
    // RegularMesh<DIM> &mesh_v = irregular_mesh_v->regularMesh();
    // RegularMesh<DIM> &mesh_p = irregular_mesh_p->regularMesh();

    // for (int i = 0; i < mesh_p.n_geometry(0); ++i)
    // {
    // 	(*mesh_p.h_geometry<0>(i))[0] = point(i)[0];
    // 	(*mesh_p.h_geometry<0>(i))[1] = point(i)[1];
    // }

    // /// 更新p网格的中点. 学习这个过程, 似乎可以去掉单元对应? 
    // for (int j = 0; j < mesh_p.n_geometry(1); ++j)
    // {
    // 	GeometryBM &bnd = mesh_p.geometry(1, j);
    // 	(*mesh_p.h_geometry<1>(bnd.index())->child[1]->vertex[0])[0]
    // 	    = 0.5 * ((*mesh_p.h_geometry<0>(bnd.vertex(0)))[0] + 
    // 		     (*mesh_p.h_geometry<0>(bnd.vertex(1)))[0]);
    // 	(*mesh_p.h_geometry<1>(bnd.index())->child[1]->vertex[0])[1]
    // 	    = 0.5 * ((*mesh_p.h_geometry<0>(bnd.vertex(0)))[1] + 
    // 		     (*mesh_p.h_geometry<0>(bnd.vertex(1)))[1]);
    // }
    // irregular_mesh_p->semiregularize();
    // irregular_mesh_p->regularize(false);
    // irregular_mesh_v->semiregularize();
    // irregular_mesh_v->regularize(false);
};

void ISOP2P1::getMonitor()
{
    RegularMesh<DIM> &regular_mesh = irregular_mesh_p->regularMesh();
    
    FEMFunction<double, DIM> u_h_p(fem_space_p);
    FEMFunction<double, DIM> v_h_p(fem_space_p);
    /// 将速度插值到压力单元上.
    Operator::L2Interpolate(v_h[0], u_h_p);
    Operator::L2Interpolate(v_h[1], v_h_p);

    FEMSpace<double, DIM>::ElementIterator the_element = fem_space_p.beginElement();
    FEMSpace<double, DIM>::ElementIterator end_element = fem_space_p.endElement();
    for (int i = 0; the_element != end_element; ++the_element) 
    {
    	int idx_p = the_element->index();
	/// 几何信息.
	double volume = the_element->templateElement().volume();
	const QuadratureInfo<DIM>& quad_info = the_element->findQuadratureInfo(2);
	std::vector<double> jacobian = the_element->local_to_global_jacobian(quad_info.quadraturePoint());
	int n_quadrature_point = quad_info.n_quadraturePoint();
	std::vector<Point<DIM> > q_point = the_element->local_to_global(quad_info.quadraturePoint());
	std::vector<std::vector<double> > vx_gradient = u_h_p.gradient(q_point, *the_element);
	std::vector<std::vector<double> > vy_gradient = v_h_p.gradient(q_point, *the_element);
	std::vector<double> vx_value = u_h_p.value(q_point, *the_element);
	std::vector<double> vy_value = v_h_p.value(q_point, *the_element);
	
	double d = 0, area = 0, vh_L2norm = 0.0;
	for (int l = 0; l < n_quadrature_point; ++l) 
	{
	    double Jxw = quad_info.weight(l) * jacobian[l] * volume;
	    area += Jxw;
	    vh_L2norm += vx_value[l] * vx_value[l] + vy_value[l] * vy_value[l] + eps; 
	    d += Jxw * (innerProduct(vx_gradient[l], vx_gradient[l]) + innerProduct(vy_gradient[l], vy_gradient[l]));
	}
	monitor(i) = d / area;
    }

    double max_monitor = *std::max_element(monitor().begin(), monitor().end());
    std::cerr << "max monitor=" << *std::max_element(monitor().begin(), monitor().end())
	      << "\tmin monitor=" << *std::min_element(monitor().begin(), monitor().end())
	      << std::endl;
    smoothMonitor(5);
    for (int i = 0; i < n_geometry(DIM); ++i)
    {
	monitor(i) = 1. / sqrt(1. + alpha * pow((monitor(i)), beta));
    }
};
void ISOP2P1::updateSolution()
{
    /// 在MovingMesh2D的moveMesh()中,实际上是先做的updateSolution,然后再更新网格,
    /// 压力网格移动步长.
    double msl = moveStepLength();
    /// 网格移动后,更新速度和压力，按照邸亚娜老师文章上的方法,
    /// 保证在update的过程中保证质量守恒条件. 实际上是求解
    /// 类似时间发展的对流方程.
    
    /// 先更新插值点.
    fem_space_p.updateDofInterpPoint();
    fem_space_v.updateDofInterpPoint();

    /// 因为网格移动了, 插值点变化, 矩阵需要重新构造.
    /// 是不是可以不用重新更改稀疏矩阵结构?
    // buildMatrixStruct();
    buildMatrix();
    
    int n_dof_v = fem_space_v.n_dof();
    int n_dof_p = fem_space_p.n_dof();
    int n_total_dof = 2 * n_dof_v + n_dof_p;

    /// 备份数值解.
    FEMFunction<double, DIM> _u_h(v_h[0]);
    FEMFunction<double, DIM> _v_h(v_h[1]);
    FEMFunction<double, DIM> _p_h(p_h);
    
    /// 虚拟时间, dt 可以取到1.
    double _dt = 1.0;
    /// 为了稳定性考虑, m至少取3, 少于3的话, 数值解会不稳定.
    /// 三阶显示Rugue-Kutta方法.
    for (int m = 3; m > 0; --m)
    {
	/// 系数矩阵直接使用 Stokes 矩阵结构.
	matrix.reinit(sp_stokes);
	/// (0, 0) 
	for (int i = 0; i < sp_vxvx.n_nonzero_elements(); ++i)
	    matrix.global_entry(index_vxvx[i]) = mat_v_mass.global_entry(i);
	/// (1, 1) 这两个对角块对应质量算子.
	for (int i = 0; i < sp_vyvy.n_nonzero_elements(); ++i)
	    matrix.global_entry(index_vyvy[i]) = mat_v_mass.global_entry(i);
	/// (0, 2) 这个不是方阵. 在矩阵结构定义的时候已经直接排除了对角元优
	/// 先.
	for (int i = 0; i < sp_pvx.n_nonzero_elements(); ++i)
	    matrix.global_entry(index_pvx[i]) = (_dt / m) * mat_pvx_divT.global_entry(i);

	/// (1, 2)
	for (int i = 0; i < sp_pvy.n_nonzero_elements(); ++i)
	    matrix.global_entry(index_pvy[i]) = (_dt / m) * mat_pvy_divT.global_entry(i);

	/// (2, 0)
	for (int i = 0; i < sp_vxp.n_nonzero_elements(); ++i)
	    matrix.global_entry(index_vxp[i]) = (_dt / m) * mat_vxp_div.global_entry(i);
	
	/// (2, 1) 这四块直接复制散度矩阵. 
	for (int i = 0; i < sp_vyp.n_nonzero_elements(); ++i)
	    matrix.global_entry(index_vyp[i]) = (_dt / m) * mat_vyp_div.global_entry(i);
	/// 问题右端项.
	rhs.reinit(n_total_dof);

	/// debug
	const std::size_t * rowstart = sp_stokes.get_rowstart_indices();
	const unsigned int * colnum = sp_stokes.get_column_numbers();
	std::ofstream mat_deb;
	mat_deb.open("mat.m", std::ofstream::out);
	mat_deb.setf(std::ios::fixed);
	mat_deb.precision(20);
	    
	mat_deb << "A = sparse(" << n_total_dof << ", " << n_total_dof << ");" << std::endl;
	for (int i = 0; i < n_total_dof; ++i)
	{
	    for (int j = rowstart[i]; j < rowstart[i + 1]; ++j)
	    {
		mat_deb << "A(" << i + 1 << ", " << colnum[j] + 1 << ")=" 
			<< matrix.global_entry(j) << ";" << std::endl;
	    }
	}
	mat_deb.close();
	std::cout << "mat output" << std::endl;
	getchar();
	/// debug

	FEMSpace<double, DIM>::ElementIterator the_element_v = fem_space_v.beginElement();
	FEMSpace<double, DIM>::ElementIterator end_element_v = fem_space_v.endElement();
	/// 遍历速度单元, 拼装相关系数矩阵和右端项.
	for (the_element_v = fem_space_v.beginElement(); 
	     the_element_v != end_element_v; ++the_element_v) 
	{
	    /// 当前单元信息.
	    double volume = the_element_v->templateElement().volume();
	    /// 积分精度, u 和 p 都是 1 次, 梯度和散度 u 都是常数. 因此矩阵拼
	    /// 装时积分精度不用超过 1 次.  
	    const QuadratureInfo<DIM>& quad_info = the_element_v->findQuadratureInfo(1);
	    std::vector<double> jacobian = the_element_v->local_to_global_jacobian(quad_info.quadraturePoint());
	    int n_quadrature_point = quad_info.n_quadraturePoint();
	    std::vector<Point<DIM> > q_point = the_element_v->local_to_global(quad_info.quadraturePoint());
	    /// 速度单元信息.
	    std::vector<std::vector<double> > basis_value_v = the_element_v->basis_function_value(q_point);
	    std::vector<std::vector<std::vector<double> > > basis_gradient_v = the_element_v->basis_function_gradient(q_point);
	    std::vector<double> vx_value = _u_h.value(q_point, *the_element_v);
	    std::vector<double> vy_value = _v_h.value(q_point, *the_element_v);
	    std::vector<double> fx_value = source_v[0].value(q_point, *the_element_v);
	    std::vector<double> fy_value = source_v[1].value(q_point, *the_element_v);
	    std::vector<std::vector<double> > vx_gradient = v_h[0].gradient(q_point, *the_element_v);
            std::vector<std::vector<double> > vy_gradient = v_h[1].gradient(q_point, *the_element_v);
	    const std::vector<int>& element_dof_v = the_element_v->dof();
	    int n_element_dof_v = the_element_v->n_dof();
	    /// 对应压力单元信息,可能用不到，先放在这里.
	    Element<double, DIM> &p_element = fem_space_p.element(index_v2p[the_element_v->index()]);
	    const std::vector<int>& element_dof_p = p_element.dof();
	    std::vector<std::vector<std::vector<double> > > basis_gradient_p = p_element.basis_function_gradient(q_point);
	    std::vector<std::vector<double> >  basis_value_p = p_element.basis_function_value(q_point);
	    int n_element_dof_p = p_element.n_dof();
	    std::vector<double> p_value = _p_h.value(q_point, p_element);
	    /// 注意这里的moveDirection实际上是速度网格上积分点的移动, 因为速度单元是在压力单元的一部分, 位置有重合.
	    /// 因此可以直接掉用 MovingMesh2D中内置计算点移动方向的函数moveDirection().
	    std::vector<std::vector<double> > move_vector = moveDirection(q_point, p_element.index());

	    for (int l = 0; l < n_quadrature_point; ++l)
	    {
		double Jxw = quad_info.weight(l) * jacobian[l] * volume;
		for (int j = 0; j < n_element_dof_v; ++j)
		{
		    double rhs_cont = Jxw * basis_value_v[j][l] * (vx_value[l]); // + (_dt / m) * msl * innerProduct(move_vector[l], vx_gradient[l]));
		    rhs(element_dof_v[j]) += rhs_cont;

		    rhs_cont = Jxw * basis_value_v[j][l] * (vy_value[l]); // + (_dt / m) * msl * innerProduct(move_vector[l], vy_gradient[l]));			
		    rhs(n_dof_v + element_dof_v[j]) += rhs_cont;
		}
	    }
	}
	    
	/// 构建系数矩阵和右端项.
	Vector<double> x(n_total_dof);
	for (int i = 0; i < n_dof_v; ++i)
	{
	    x(i) = v_h[0](i);
	    x(i + n_dof_v) = v_h[1](i);
	}
	
	for (int i = 0; i < n_dof_p; ++i)
	    x(i + 2 * n_dof_v) = p_h(i);
	/// 处理边界.
	boundaryValueUpdateSolution(x);

	clock_t t_cost = clock();
	dealii::SolverControl solver_control (4000000, l_Euler_tol, check);
	/// AMG 预处理不行,试一下ILU预处理.
	/// 不完全LU分解.     
	dealii::SparseILU <double> preconditioner;
	preconditioner.initialize(matrix);

	SolverMinRes<Vector<double> > minres(solver_control);
	minres.solve (matrix, x, rhs, PreconditionIdentity());
	t_cost = clock() - t_cost;

	std::cout << "time cost: " << (((float)t_cost) / CLOCKS_PER_SEC) << std::endl;
	///速度数值解和压力数值解.
	for (int i = 0; i < n_dof_v; ++i)
	{
	    v_h[0](i) = x(i);
	    v_h[1](i) = x(i + n_dof_v);
	}	
	for (int i = 0; i < n_dof_p; ++i)
	    p_h(i) = x(i + 2 * n_dof_v);
    }
    outputTecplot("NS_Euler");
    outputSolution();
}; 

void ISOP2P1::updateMesh()
{
    for (int i = 0; i < n_geometry(0); ++i) 
    {
	point(i)[0] += moveStepLength() * moveDirection(i)[0];
	point(i)[1] += moveStepLength() * moveDirection(i)[1];
    }
    /// 移动irregular_mesh_p树上的网格点.
    RegularMesh<DIM> &mesh_p = irregular_mesh_p->regularMesh();
    for (int i = 0; i < mesh_p.n_geometry(0); ++i)
    {
    	(*mesh_p.h_geometry<0>(i))[0] = point(i)[0];
    	(*mesh_p.h_geometry<0>(i))[1] = point(i)[1];
    }
    /// 更新p网格的中点.
    for (int j = 0; j < mesh_p.n_geometry(1); ++j)
    {
    	GeometryBM &bnd = mesh_p.geometry(1, j);
    	(*mesh_p.h_geometry<1>(bnd.index())->child[1]->vertex[0])[0]
    	    = 0.5 * ((*mesh_p.h_geometry<0>(bnd.vertex(0)))[0] + 
    		     (*mesh_p.h_geometry<0>(bnd.vertex(1)))[0]);
    	(*mesh_p.h_geometry<1>(bnd.index())->child[1]->vertex[0])[1]
    	    = 0.5 * ((*mesh_p.h_geometry<0>(bnd.vertex(0)))[1] + 
    		     (*mesh_p.h_geometry<0>(bnd.vertex(1)))[1]);
    }

    // /// irregular_mesh_p 和 irregular_mesh_v 在同一棵树上, 因此, 移动了 irregular_mesh_v, 
    /// irregular_mesh_p上的网格点也移动了, 不能重复移动两次.
    irregular_mesh_p->semiregularize();
    irregular_mesh_p->regularize(false);
    irregular_mesh_v->semiregularize();
    irregular_mesh_v->regularize(false);
    /// 更新一下插值点.
    fem_space_v.updateDofInterpPoint();
    fem_space_p.updateDofInterpPoint();
}
void ISOP2P1::outputSolution()
{
    p_h.writeOpenDXData("p_h.dx");
    outputPhysicalMesh("E");
};

void ISOP2P1::movingMesh()
{
    // /// 修改最大移动次数.
    // maxStep() = max_step;
    // /// 设置容忍量.
    // tolerence() = 0.01;

    // moveMesh();
    // /// 在一步网格移动之后, 同步一下速度和压力网格.
    // syncMesh();
    // fem_space_p.updateDofInterpPoint();
    // fem_space_v.updateDofInterpPoint();
};