From: Shengli Xu <shengli.xu@gm...>  20080322 11:54:29

Hi Ben and Roy In libMesh 0.5.0, The assemble function runs very quickly. But in libMesh0.6.2, It is very slowly. The memory used in libMesh0.5.0 is much larger than in libMesh0.6.2 before call assemble function. The program is the same. Why? I debugged my programe. The function of apply_periodic_bc() and add_system_matrix() are both called in es.init(), but apply_periodic_bc() is in the front of add_system_matrix(). I think that the function "this>init_matrix()" in void ImplicitSystem::init_data() should consider the information of periodic b.c. which is in DofMap. The program structure is like this: int main() { EquationSystems es(mesh); LinearImplicitSystem &system = ... .... system.add_variable("u", FIRST); system.add_variable("v", FIRST); system.add_variable("w",FIRST); system.attach_assemble_function(assemble_elastic); system.attach_constraint_function(apply_periodic_bc_elastic); es.init(); .... system.solve(); ... } void assemble_elastic(EquationSystems& es, const std::string& system_name) { ... for(; el!=end_el; ++el) { // Gauss integrate { .... } dof_map.constrain_element_matrix_and_vector(Ke,Fe,dof_indices); system.matrix>add_matrix(Ke,dof_indices); system.rhs>add_vector(Fe,dof_indices); } // end el loop } /* * Apply periodic b.c. for elastic homogenization * Read the node couple to apply periodic b.c. */ void apply_periodic_bc_elastic(EquationSystems& es, const std::string& system_name) { std::cout<<"Apply periodic boundary condition for elastic homogenization..."<<std::endl; Mesh &mesh=es.get_mesh(); LinearImplicitSystem& system=es.get_system <LinearImplicitSystem>(system_name.c_str()); DofMap& dof_map = system.get_dof_map(); std::cout<<"Imput file of periodic b.c. node for elastic module"<<std::endl; std::string PeribcFile; std::cin>>PeribcFile; // The file of node couple std::cout<<"PeribcFile="<<PeribcFile<<std::endl; std::ifstream pbcfile(PeribcFile.c_str()); unsigned int p2num, p4num; std::vector<int> p2_1, p2_2, p4_1, p4_2, p4_3, p4_4; pbcfile>>p2num; for(unsigned int i=0; i<p2num; i++) { int m,n; pbcfile>>m>>n; // Read the node number p2_1.push_back(m1); p2_2.push_back(n1); } pbcfile>>p4num; for(unsigned int i=0; i<p4num; i++) { int m1, m2, m3, m4; pbcfile>>m1>>m2>>m3>>m4; p4_1.push_back(m11); p4_2.push_back(m21); p4_3.push_back(m31); p4_4.push_back(m41); } pbcfile.close(); // for(unsigned int i=0;i<p2num;i++){ unsigned int n1=p2_1.at(i); unsigned int n2=p2_2.at(i); const Node &curr_n1 = mesh.node(n1); const Node &curr_n2 = mesh.node(n2); unsigned int u1 = curr_n1.dof_number(0,0,0); unsigned int v1 = curr_n1.dof_number(0,1,0); unsigned int w1 = curr_n1.dof_number(0,2,0); unsigned int u2 = curr_n2.dof_number(0,0,0); unsigned int v2 = curr_n2.dof_number(0,1,0); unsigned int w2 = curr_n2.dof_number(0,2,0); DofConstraintRow cu,cv,cw; cu[u1]=1.0; cv[v1]=1.0; cw[w1]=1.0; dof_map.add_constraint_row(u2,cu); dof_map.add_constraint_row(v2,cv); dof_map.add_constraint_row(w2,cw); } // for(unsigned int i=0;i<p4num;i++){ unsigned int n1=p4_1.at(i); unsigned int n2=p4_2.at(i); unsigned int n3=p4_3.at(i); unsigned int n4=p4_4.at(i); const Node &curr_n1=mesh.node(n1); const Node &curr_n2=mesh.node(n2); const Node &curr_n3=mesh.node(n3); const Node &curr_n4=mesh.node(n4); unsigned int u1=curr_n1.dof_number(0,0,0); unsigned int v1=curr_n1.dof_number(0,1,0); unsigned int w1=curr_n1.dof_number(0,2,0); unsigned int u2=curr_n2.dof_number(0,0,0); unsigned int v2=curr_n2.dof_number(0,1,0); unsigned int w2=curr_n2.dof_number(0,2,0); unsigned int u3=curr_n3.dof_number(0,0,0); unsigned int v3=curr_n3.dof_number(0,1,0); unsigned int w3=curr_n3.dof_number(0,2,0); unsigned int u4=curr_n4.dof_number(0,0,0); unsigned int v4=curr_n4.dof_number(0,1,0); unsigned int w4=curr_n4.dof_number(0,2,0); DofConstraintRow cu,cv,cw; cu[u1]=1.0; cv[v1]=1.0; cw[w1]=1.0; dof_map.add_constraint_row(u2,cu); dof_map.add_constraint_row(u3,cu); dof_map.add_constraint_row(u4,cu); dof_map.add_constraint_row(v2,cv); dof_map.add_constraint_row(v3,cv); dof_map.add_constraint_row(v4,cv); dof_map.add_constraint_row(w2,cw); dof_map.add_constraint_row(w3,cw); dof_map.add_constraint_row(w4,cw); } } Sincerely Shengli 2008/3/20, Benjamin Kirk <benjamin.kirk@...>: > > > You must be adding the constraint rows after the matrix has been > > preallocated... Regardless of whether PeriodicBoundary calls would > > fix your problem, we ought to make sure that users are able and > > encouraged to do efficient preallocation for arbitrary user > > constraints, too. > > The DOF connectivity is established during EquationSystems::init(). So > you must add any userdefined constraints before this time. The matrix > itself is not actually allocated until the first call to > EquationSystems::solve(), but when it is allocated it is using a > precomputed sparsity pattern. > > The proper place to do this is to provide a function pointer which is > called by System::user_constrain(): > > http://libmesh.sourceforge.net/doxygen/classSystem.php#8dd9e548cdc41149155da0a92eefa591 > > Try that, it should fix your problem. > > Ben >  Best regards, Yours sincerely ShengliXu Department of Engineering Mechanics State Key Laboratory of Structural Analysis for Industrial Equipment Dalian University of Technology Dalian, 116023, P. R. China Email: shengli.xu.xu@... ========================== 