From: Omar Al-A. <oma...@gm...> - 2011-05-10 15:33:17
|
Hi everyone, I am trying to define the boundary conditions in a way which I am more familiar with. I have only one type of BC which is Drichilet and I am doing that simply by making indices that are on the outer surface of the domain to be zero in the stiffness (k) and mass (m) matrix and 1 for the mass matrix if index i==j. Ok I am trying to do that, but I am facing a problem for the mass matrix when the index i==j, since I am using: matrix_B.add_matrix (Me, dof_indices); which will add several 1 to the diagonal terms in the matrix. I know that I am applying my BC in a very poor way, I was trying to imitate one of the example to my purpose and this what I am doing: for (unsigned int side=0; side<elem->n_sides(); side++){ if (elem->neighbor(side) == NULL) { const std::vector<std::vector<Real> >& phi_face = fe_face->get_phi(); const std::vector<Real>& JxW_face = fe_face->get_JxW(); const std::vector<Point >& qface_point = fe_face->get_xyz(); // Compute the shape function values on the element // face. fe_face->reinit(elem, side); // Loop over the face quadrature points for integration. // for (unsigned int qp=0; qp<qface.n_points(); qp++) // { // Matrix contribution of the L2 projection. for (unsigned int i=0; i<phi_face.size(); i++){ for (unsigned int j=0; j<phi_face.size(); j++){ Ke(i,j) = 0.; if (i==j){Me(i,j) = 1.;} else {Me(i,j) = 0.;} } // } // end face quadrature point loop } // end if (elem->neighbor(side) == NULL) } // end boundary condition section } matrix_A.add_matrix (Ke, dof_indices); matrix_B.add_matrix (Me, dof_indices); Appreciate if someone can find me a way the will make the diagonal terms of the surface of mass matrix to be equal 1. And if someone has better suggestion than the miss that I am doing in ( for (unsigned int i=0; i<phi_face.size(); i++){ for (unsigned int j=0; j<phi_face.size(); j++){ ......) I will appreciate that too. Thanks in advance, Omar |
From: Boyce G. <gri...@ci...> - 2011-05-10 19:01:03
|
On 5/10/11 11:33 AM, Omar Al-Abbasi wrote: > Hi everyone, > > I am trying to define the boundary conditions in a way which I am more > familiar with. I have only one type of BC which is Drichilet and I am doing > that simply by making indices that are on the outer surface of the domain to > be zero in the stiffness (k) and mass (m) matrix and 1 for the mass matrix > if index i==j. > > Ok I am trying to do that, but I am facing a problem for the mass matrix > when the index i==j, since I am using: > matrix_B.add_matrix (Me, dof_indices); > > which will add several 1 to the diagonal terms in the matrix. I think that the easiest way to do this is to make two sweeps through the mesh. In the first sweep, set the matrix elements without worrying about boundary conditions. Once you have finished this initial assembly, call SparseMatrix::close() to "flush" the matrix coefficients. Then make a second sweep through the mesh and reset the appropriate entries by calling SparseMatrix::set(). -- Boyce |
From: John P. <jwp...@gm...> - 2011-05-10 19:04:44
|
On Tue, May 10, 2011 at 1:00 PM, Boyce Griffith <gri...@ci...> wrote: > > > On 5/10/11 11:33 AM, Omar Al-Abbasi wrote: >> Hi everyone, >> >> I am trying to define the boundary conditions in a way which I am more >> familiar with. I have only one type of BC which is Drichilet and I am doing >> that simply by making indices that are on the outer surface of the domain to >> be zero in the stiffness (k) and mass (m) matrix and 1 for the mass matrix >> if index i==j. >> >> Ok I am trying to do that, but I am facing a problem for the mass matrix >> when the index i==j, since I am using: >> matrix_B.add_matrix (Me, dof_indices); >> >> which will add several 1 to the diagonal terms in the matrix. > > I think that the easiest way to do this is to make two sweeps through > the mesh. In the first sweep, set the matrix elements without worrying > about boundary conditions. Once you have finished this initial > assembly, call SparseMatrix::close() to "flush" the matrix coefficients. > Then make a second sweep through the mesh and reset the appropriate > entries by calling SparseMatrix::set(). I think DenseMatrix::condense() is the appropriate way to handle this. We don't have an example which does this, but the interface is relatively straightforward, and is defined in include/numerics/dense_matrix_base.h. -- John |
From: Boyce G. <gri...@ci...> - 2011-05-10 19:21:19
|
On 5/10/11 3:04 PM, John Peterson wrote: > On Tue, May 10, 2011 at 1:00 PM, Boyce Griffith<gri...@ci...> wrote: >> >> >> On 5/10/11 11:33 AM, Omar Al-Abbasi wrote: >>> Hi everyone, >>> >>> I am trying to define the boundary conditions in a way which I am more >>> familiar with. I have only one type of BC which is Drichilet and I am doing >>> that simply by making indices that are on the outer surface of the domain to >>> be zero in the stiffness (k) and mass (m) matrix and 1 for the mass matrix >>> if index i==j. >>> >>> Ok I am trying to do that, but I am facing a problem for the mass matrix >>> when the index i==j, since I am using: >>> matrix_B.add_matrix (Me, dof_indices); >>> >>> which will add several 1 to the diagonal terms in the matrix. >> >> I think that the easiest way to do this is to make two sweeps through >> the mesh. In the first sweep, set the matrix elements without worrying >> about boundary conditions. Once you have finished this initial >> assembly, call SparseMatrix::close() to "flush" the matrix coefficients. >> Then make a second sweep through the mesh and reset the appropriate >> entries by calling SparseMatrix::set(). > > I think DenseMatrix::condense() is the appropriate way to handle this. > > We don't have an example which does this, but the interface is > relatively straightforward, and is defined in > include/numerics/dense_matrix_base.h. I think DenseMatrix::condense() will work to zero out values, but if you want to put some nonzero value in, e.g., the diagonal of the assembled matrix, it seems like you need to do that after an initial assembly, don't you? -- Boyce |
From: Kirk, B. (JSC-EG311) <ben...@na...> - 2011-05-10 19:26:41
|
> I think DenseMatrix::condense() will work to zero out values, but if you > want to put some nonzero value in, e.g., the diagonal of the assembled > matrix, it seems like you need to do that after an initial assembly, > don't you? Not sure if this is the same concern as the previous - about adding the diagonal multiple times, but if it is that is not a problem because while you do add the diagonal multiple times you also add the RHS multiple times. So for a DOF shared by say 5 elements rather than [0 ... 0 1 0 ....0] dof = val you get instead [0 ... 0 5 0 ....0] dof = 5*val -Ben |
From: Omar Al-A. <oma...@gm...> - 2011-05-10 23:10:45
|
Thank you all for all the suggestions you made. In fact I tried doing two seperate sweeps as Boyce suggested, but it seems that I am doing a mistake which leads to no change in mass matrix. This how I am doing that after the first sweep for assembling the matrices: . . . for ( ; el != end_el; ++el) { const Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); fe->reinit (elem); Ke.resize (dof_indices.size(), dof_indices.size()); Me.resize (dof_indices.size(), dof_indices.size()); for (unsigned int side=0; side<elem->n_sides(); side++){ if (elem->neighbor(side) == NULL) { const std::vector<std::vector<Real> >& phi_face = fe_face->get_phi(); fe_face->reinit(elem, side); for (unsigned int i=0; i<phi_face.size(); i++){ matrix_B.set(dof_indices[i],dof_indices[i],1.); } // end if (elem->neighbor(side) == NULL) } // end boundary condition section } } Also, what is the best way to loop over the nodes which live on the outer surface of the domain? Is it just the way I am doing it for (unsigned int side=0; side<elem->n_sides(); side++){ if (elem->neighbor(side) == NULL) and then for (unsigned int i=0; i<phi_face.size(); i++){ .... Thanks again, Omar |
From: Derek G. <fri...@gm...> - 2011-05-11 14:54:02
|
The easiest way to do what you want is by using matrix.zero_rows(). The way it works is that you fill up a vector of dof indices that correspond to rows that you want to zero and put a 1 on the diagonal. Then you just call matrix.zero_rows(dofs_to_zero, 1). The second argument to the function is the number you want to place on the diagonal (zero by default). So.... fill up your matrix like normal, keeping track of rows you want to wipe out. Close the matrix. Call zero_rows(). Close the matrix again. A quick note on this... if you are using Petsc this will modify your nonzero pattern (ie remove all off-diagonal nonzero entries for these rows). That's not normally what you want because the next time you come around to fill the matrix and you do some insertions on that row it will be incredibly slow (because Petsc will have to reallocate those entries). You can get beyond this by using something like: MatSetOption(static_cast<PetscMatrix<Number> &>(my_matrix).mat(), MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); That is the correct line for Petsc 3.1. If you are using a different version there will most likely be a slightly different way to tell it to keep the nonzero pattern. I recommend setting this each time before you start filling the matrix up. Petsc seems to like to forget about this option..... Derek On May 10, 2011, at 5:10 PM, Omar Al-Abbasi wrote: > Thank you all for all the suggestions you made. > In fact I tried doing two seperate sweeps as Boyce suggested, but it seems > that I am doing a mistake which leads to no change in mass matrix. > > This how I am doing that after the first sweep for assembling the matrices: > . > . > . > for ( ; el != end_el; ++el) > { > const Elem* elem = *el; > dof_map.dof_indices (elem, dof_indices); > fe->reinit (elem); > Ke.resize (dof_indices.size(), dof_indices.size()); > Me.resize (dof_indices.size(), dof_indices.size()); > > for (unsigned int side=0; side<elem->n_sides(); side++){ > if (elem->neighbor(side) == NULL) > { > const std::vector<std::vector<Real> >& phi_face = > fe_face->get_phi(); > fe_face->reinit(elem, side); > for (unsigned int i=0; i<phi_face.size(); i++){ > matrix_B.set(dof_indices[i],dof_indices[i],1.); > } // end if (elem->neighbor(side) == NULL) > } // end boundary condition section > } > } > > Also, what is the best way to loop over the nodes which live on the outer > surface of the domain? > Is it just the way I am doing it for (unsigned int side=0; > side<elem->n_sides(); side++){ > if (elem->neighbor(side) == NULL) > and then for (unsigned int i=0; i<phi_face.size(); i++){ .... > > Thanks again, > Omar > ------------------------------------------------------------------------------ > Achieve unprecedented app performance and reliability > What every C/C++ and Fortran developer should know. > Learn how Intel has extended the reach of its next-generation tools > to help boost performance applications - inlcuding clusters. > http://p.sf.net/sfu/intel-dev2devmay > _______________________________________________ > Libmesh-users mailing list > Lib...@li... > https://lists.sourceforge.net/lists/listinfo/libmesh-users |
From: Omar Al-A. <oma...@gm...> - 2011-05-11 17:31:25
|
Brilliant, that would do the job in a very easy way. Ok last question that I would ask related to this topic :-) I not sure how to find the nodes (DOFs) which live on the boundary (outer surface of my 3D domain). I appreciate if someone can finalize this issue for me. Kindest regards, Omar On Wed, May 11, 2011 at 10:53 AM, Derek Gaston <fri...@gm...> wrote: > The easiest way to do what you want is by using matrix.zero_rows(). The > way it works is that you fill up a vector of dof indices that correspond to > rows that you want to zero and put a 1 on the diagonal. Then you just call > matrix.zero_rows(dofs_to_zero, 1). The second argument to the function is > the number you want to place on the diagonal (zero by default). > > So.... fill up your matrix like normal, keeping track of rows you want to > wipe out. Close the matrix. Call zero_rows(). Close the matrix again. > > A quick note on this... if you are using Petsc this will modify your > nonzero pattern (ie remove all off-diagonal nonzero entries for these rows). > That's not normally what you want because the next time you come around to > fill the matrix and you do some insertions on that row it will be incredibly > slow (because Petsc will have to reallocate those entries). You can get > beyond this by using something like: > > MatSetOption(static_cast<PetscMatrix<Number> &>(my_matrix).mat(), > MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); > > That is the correct line for Petsc 3.1. If you are using a different > version there will most likely be a slightly different way to tell it to > keep the nonzero pattern. I recommend setting this each time before you > start filling the matrix up. Petsc seems to like to forget about this > option..... > > Derek > > On May 10, 2011, at 5:10 PM, Omar Al-Abbasi wrote: > > > Thank you all for all the suggestions you made. > > In fact I tried doing two seperate sweeps as Boyce suggested, but it > seems > > that I am doing a mistake which leads to no change in mass matrix. > > > > This how I am doing that after the first sweep for assembling the > matrices: > > . > > . > > . > > for ( ; el != end_el; ++el) > > { > > const Elem* elem = *el; > > dof_map.dof_indices (elem, dof_indices); > > fe->reinit (elem); > > Ke.resize (dof_indices.size(), dof_indices.size()); > > Me.resize (dof_indices.size(), dof_indices.size()); > > > > for (unsigned int side=0; side<elem->n_sides(); side++){ > > if (elem->neighbor(side) == NULL) > > { > > const std::vector<std::vector<Real> >& phi_face = > > fe_face->get_phi(); > > fe_face->reinit(elem, side); > > for (unsigned int i=0; i<phi_face.size(); i++){ > > matrix_B.set(dof_indices[i],dof_indices[i],1.); > > } // end if (elem->neighbor(side) == NULL) > > } // end boundary condition section > > } > > } > > > > Also, what is the best way to loop over the nodes which live on the outer > > surface of the domain? > > Is it just the way I am doing it for (unsigned int side=0; > > side<elem->n_sides(); side++){ > > if (elem->neighbor(side) == NULL) > > and then for (unsigned int i=0; i<phi_face.size(); i++){ .... > > > > Thanks again, > > Omar > > > ------------------------------------------------------------------------------ > > Achieve unprecedented app performance and reliability > > What every C/C++ and Fortran developer should know. > > Learn how Intel has extended the reach of its next-generation tools > > to help boost performance applications - inlcuding clusters. > > http://p.sf.net/sfu/intel-dev2devmay > > _______________________________________________ > > Libmesh-users mailing list > > Lib...@li... > > https://lists.sourceforge.net/lists/listinfo/libmesh-users > > |