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(m-1);
p2_2.push_back(n-1);
}
pbcfile>>p4num;
for(unsigned int i=0; i<p4num; i++)
{
int m1, m2, m3, m4;
pbcfile>>m1>>m2>>m3>>m4;
p4_1.push_back(m1-1);
p4_2.push_back(m2-1);
p4_3.push_back(m3-1);
p4_4.push_back(m4-1);
}
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 <[EMAIL PROTECTED]>:
>
> > 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 user-defined 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
> pre-computed 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: [EMAIL PROTECTED]
==========================
-------------------------------------------------------------------------
This SF.net email is sponsored by: Microsoft
Defy all challenges. Microsoft(R) Visual Studio 2008.
http://clk.atdmt.com/MRT/go/vse0120000070mrt/direct/01/
_______________________________________________
Libmesh-users mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/libmesh-users
