[deal.II] How to use sacado package

[Jaekwang Kim]

I'd like to use sacado package for automatic differentiation. 


However, I have installed deal.ii on my computer with trilinos 
configuration set off. 


I looked back deal.ii installation manual and got to know that I have to 
install trilnos package first. 


However, when I try to install trilnos following instruction, 


cd trilinos-12.4.2
mkdir build
cd build

cmake\
-DTrilinos_ENABLE_Amesos=ON  \
-DTrilinos_ENABLE_Epetra=ON  \
-DTrilinos_ENABLE_Ifpack=ON  \
-DTrilinos_ENABLE_AztecOO=ON \
-DTrilinos_ENABLE_Sacado=ON  \
-DTrilinos_ENABLE_Teuchos=ON \
-DTrilinos_ENABLE_MueLu=ON   \
-DTrilinos_ENABLE_ML=ON  \
-DTrilinos_VERBOSE_CONFIGURE=OFF \
-DTPL_ENABLE_MPI=ON  \
-DBUILD_SHARED_LIBS=ON   \
-DCMAKE_VERBOSE_MAKEFILE=OFF \
-DCMAKE_BUILD_TYPE=RELEASE   \
-DCMAKE_INSTALL_PREFIX:PATH=$HOME/share/trilinos \
../

make install


I got error message as follows...




"

The Fortran compiler identification is unknown

CMake Error at 
cmake/tribits/core/package_arch/TribitsGlobalMacros.cmake:1638 
(ENABLE_LANGUAGE):

  No CMAKE_Fortran_COMPILER could be found.


  Tell CMake where to find the compiler by setting either the environment

  variable "FC" or the CMake cache entry CMAKE_Fortran_COMPILER to the full

  path to the compiler, or to the compiler name if it is in the PATH.

Call Stack (most recent call first):

  cmake/tribits/core/package_arch/TribitsProjectImpl.cmake:188 
(TRIBITS_SETUP_ENV)

  cmake/tribits/core/package_arch/TribitsProject.cmake:93 
(TRIBITS_PROJECT_IMPL)

  CMakeLists.txt:89 (TRIBITS_PROJECT)

"


What's the problem? 

and I would very thank if someone tell me how to change deal.ii 
configuration to set on Trilinos after that


Thank you! 


-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Dear Daniel,

Please look at the pdf file attached. Thanks

On Friday, October 21, 2016 at 4:51:22 PM UTC-5, Daniel Arndt wrote:
>
> Hamed,
>
> At least I can't open your file in a way to clearly see what your 
> bilinearform is, but could just guess from the implementation. Can you 
> either write it here (using LaTeX or similar) or upload the file as pdf 
> again?
>
> Best,
> Daniel
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


Bilinearform.pdf
Description: Adobe PDF document

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Daniel Arndt]

Hamed,

At least I can't open your file in a way to clearly see what your 
bilinearform is, but could just guess from the implementation. Can you 
either write it here (using LaTeX or similar) or upload the file as pdf 
again?

Best,
Daniel

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Error in GridIn::read_abaqus for a mesh with two blocks

[Oded Yaakobi]

 

Hi Jean-Paul,

 

I wanted to verify that the issue that I encountered is not related to the 
fact that I am using a version of dealii which is not the most updated one, 
so I installed the current version of dealii from git. I run two versions 
of my code with ‘apply_all_indicators_to_manifolds’ either false or true - 
see the attached code. In either case I got similar, although not exactly 
the same output. The value of the boundary indicator that is noted 
regarding the problematic quad is 1 or 0 respectively.

 

For ‘apply_all_indicators_to_manifolds’ = false:

 

Reading block 1 file:

block 1 file was read successfully

Reading block 2 file:

block 2 file was read successfully

Reading block 1 and 2 file:

Exception on processing internal UCD data: 

 



An error occurred in line <3008> of file 

 
in function

static void 
dealii::internal::Triangulation::Implementation::create_triangulation(const 
std::vector > &, const std::vector > &, const 
dealii::SubCellData &, Triangulation<3, spacedim> &) [spacedim = 3]

The violated condition was: 

quad->at_boundary()

Additional information: 

The input data for creating a triangulation contained information about 
a quad with indices 1289, 1398, 1288, and 1397 that is described to have 
boundary indicator 1. However, this is an internal quad not located on the 
boundary. You cannot assign a boundary indicator to it.

 

If this happened at a place where you call 
Triangulation::create_triangulation() yourself, you need to check the 
SubCellData object you pass to this function.

 

If this happened in a place where you are reading a mesh from a file, then 
you need to investigate why such a quad ended up in the input file. A 
typical case is a geometry that consisted of multiple parts and for which 
the mesh generator program assumes that the interface between two parts is 
a boundary when that isn't supposed to be the case, or where the mesh 
generator simply assigns 'geometry indicators' to quads at the surface of a 
part that are not supposed to be interpreted as 'boundary indicators'.



 

 

 



Exception on processing: 

 



An error occurred in line <899> of file 

 
in function

void dealii::GridIn<3, 3>::read_abaqus(std::istream &, const bool) [dim 
= 3, spacedim = 3]

The violated condition was: 

false

Additional information: 

Internal conversion from ABAQUS file to UCD format was unsuccessful. 
More information is provided in an error message printed above. Are you 
sure that your ABAQUS mesh file conforms with the requirements listed in 
the documentation?



 

Aborting!



 

 

For ‘apply_all_indicators_to_manifolds’ = true:

 

Reading block 1 file:

block 1 file was read successfully

Reading block 2 file:

block 2 file was read successfully

Reading block 1 and 2 file:

Exception on processing internal UCD data: 

 



An error occurred in line <3008> of file 

 
in function

static void 
dealii::internal::Triangulation::Implementation::create_triangulation(const 
std::vector > &, const std::vector > &, const 
dealii::SubCellData &, Triangulation<3, spacedim> &) [spacedim = 3]

The violated condition was: 

quad->at_boundary()

Additional information: 

The input data for creating a triangulation contained information about 
a quad with indices 1289, 1398, 1288, and 1397 that is described to have 
boundary indicator 0. However, this is an internal quad not located on the 
boundary. You cannot assign a boundary indicator to it.

 

If this happened at a place where you call 
Triangulation::create_triangulation() yourself, you need to check the 
SubCellData object you pass to this function.

 

If this happened in a place where you are reading a mesh from a file, then 
you need to investigate why such a quad ended up in the input file. A 
typical case is a geometry that consisted of multiple parts and for which 
the mesh generator program assumes that the interface between two parts is 
a boundary when that isn't supposed to be the case, or where the mesh 
generator simply assigns 'geometry indicators' to quads at the surface of a 
part that are not supposed to be interpreted as 'boundary indicators'.



 

 

 



Exception on processing: 

 



An error occurred in line <899> of file 

 
in function

void dealii::GridIn<3, 3>::read_abaqus(std::istream &, const bool) [dim 
= 3, spacedim = 3]

The violated condition was: 

false

Additional infor

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Dear Wolfgang,

In attach you shall find the bilinear form and the assemble part of the 
code as well.

Thanks,

On Friday, October 21, 2016 at 11:14:52 AM UTC-5, Wolfgang Bangerth wrote:
>
>
> > On Thursday, October 20, 2016 at 9:48:55 PM UTC-5, Wolfgang Bangerth 
> wrote: 
> > 
> > What is the norm of the first right hand side? You now set the 
> tolerance in 
> > 
> > the first iteration to a fixed value of 1e-16, but how does this 
> compare to 
> > the previous value of 1e-12*system_rhs.l2_norm()? 
> > 
> > 
> > The system_rhs.l2_norm() before the first call to CGSolverat in the zero 
> time 
> > step was 1.82492e-09, so that the 1e-12*system_rhs.l2_norm() became in 
> the 
> > order of 1e-21 by which CGSOlver didn't converge. 
>
> OK. I have no suggestion, though, why one or the other would lead to 
> convergence or non-convergence of CG. 
>
>
> > That still seems wrong to me. As I mentioned in a previous email, 
> you 
> > ought to 
> > make sure that the matrix you build is the matrix you *want* to 
> build. 
> > 
> > 
> > I was wondering how to make sure the matrix I built is what I meant. 
>
> You need to compare properties of your bilinear form with those of your 
> matrix. For example, if your bilinear form is symmetric, then your matrix 
> is 
> as well -- and that is something you can test. If your bilinear form is 
> positive definite, then your matrix needs to be as well and as a 
> consequence 
> CG needs to converge. If CG does not seem to converge but your bilinear 
> form 
> is s.p.d., you know you have a bug in your matrix assembly. If CG 
> converges on 
> one processor but not on two for the same problem, then you know you have 
> a 
> bug in the parallel matrix assembly. Etc. 
>
> Since you haven't stated your bilinear form or the problem you are 
> solving, I 
> can't tell you whether the matrix is supposed to be s.p.d. 
>
> Best 
>   W. 
>
>
> -- 
>  
> Wolfgang Bangerth  email: bang...@colostate.edu 
>  
> www: http://www.math.colostate.edu/~bangerth/ 
>
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.


Bilinearform.docx
Description: MS-Word 2007 document

[deal.II] Re: Error in GridIn::read_abaqus for a mesh with two blocks

[Jean-Paul Pelteret]

Hi Oded,

I wasn't able to read in the .cub file (clearly its not backwards 
compatible). I did manage to import the .inp file and nothing looks amiss. 
Could you try to set the apply_all_indicators_to_manifolds 

 flag 
to true? Its possible that one of the boundary id's is being set to an 
internal face. I've had this bug crop up on occasion for unstructured 
meshes, but never for structured ones. If this fixes the problem then 
you'll have to recreate all of the boundary ID's yourself for the time 
being until the bug is fixed.

Regards,
J-P


On Friday, October 21, 2016 at 6:27:07 PM UTC+2, Oded Yaakobi wrote:
>
> Hi Jean-Paul,
>
>  
>
> I changed one of the options in the Cubit 15.2 GUI, and now 
> 'gen_input_file' does not appear in the journal file. Reading in my code 
> the resulting abaqus files behaves as before, i.e. this is the output:
>
>  
>
> Reading block 1 file:
>
> block 1 file was read successfully
>
> Reading block 2 file:
>
> block 2 file was read successfully
>
> Reading block 1 and 2 file:
>
>  
>
>  
>
> 
>
> Exception on processing: 
>
>  
>
> 
>
> An error occurred in line <864> of file 
>  in 
> function
>
> void dealii::GridIn<3, 3>::read_abaqus(std::istream &) [dim = 3, 
> spacedim = 3]
>
> The violated condition was: 
>
> false
>
> The name and call sequence of the exception was:
>
> ExcMessage("Internal conversion from ABAQUS file to UCD format was 
> unsuccessful.Are you sure that your 
> ABAQUS mesh file conforms with the 
> requirementslisted in the 
> documentation?")
>
> Additional Information: 
>
> Internal conversion from ABAQUS file to UCD format was 
> unsuccessful.Are you sure that your 
> ABAQUS mesh file conforms with the 
> requirementslisted in the documentation?
>
> 
>
>  
>
> Aborting!
>
> 
>
>  
>
> For your convenience, I attach here all the files. I hope that you would 
> be able to run the new journal file on your version of Cubit.
>
>  
>
> Thank you, 
>
> Oded
>
> On Friday, October 21, 2016 at 11:09:01 AM UTC-4, Oded Yaakobi wrote:
>>
>> Hi Jean-Paul,
>>
>>  
>>
>> Attached is the CUB file of the two-block mesh that I created.
>>
>>  
>>
>> I created the mesh using only Cubit 15.2 GUI. In saving the abaqus file, 
>> I followed the instruction of GridIn::read_abaqus. The journal file that I 
>> attached before reflects the commands that I executed through the GUI. 
>> Probably 'gen_input_file' is not defined in the old version of Cubit that 
>> you have. 
>>
>>  
>>
>> According to Cubit 15.2 documentation, 'gen_input_file' is an input file 
>> with the given file name that will be generated when the Sculpt Parallel 
>> command will be executed. This is a text file containing all sculpt options 
>> used in the command. The input file is intended to be used for batch 
>> execution of sculpt. Therefore, in my case, in which I used the GUI, 
>> 'gen_input_file' is actually an output file and not an input file.
>>
>>  
>>
>> Best,
>>
>> Oded
>>
>>
>> On Friday, October 21, 2016 at 3:40:32 AM UTC-4, Jean-Paul Pelteret wrote:
>>>
>>> Hi Oded,
>>>
>>> Two things:
>>> 1. The documentation originally stated that "Files generated using Cubit 
>>> 11.x, 12.x and 13.x are valid", but I see that its now been updated to 
>>> include 14.x and 15.x. I personally don't know whether this works with 
>>> newer versions of Cubit, but hopefully it has indeed been tested.
>>> 2. You've not provided all of the data necessary to recreate your 
>>> geometry:
>>>
>>>  CUBIT> create sphere radius 1
 Successfully created sphere volume 1 
 Journaled Command: create sphere radius 1
 CUBIT> sculpt parallel volume all processors 8 fileroot 
 "sculpt_parallel" gen_input_file "sculpt_parallel.i" overwrite xintervals 
 8 
 yintervals 8 zintervals 8 box location position -3 -3 -3 location position 
 3 3 3 smooth 8 num_laplace 2 max_opt_iters 5 opt_threshold 0.6 
 max_pcol_iters 100 pcol_threshold 0.2 sideset 2 void adapt_type 1 
 adapt_levels 2 combine import show clean
 ERROR: Unrecognized Identifier: 'gen_input_file'
>>>
>>>
>>> Can you please provide all whatever files are necessary to recreate the 
>>> geometry? Either that or provide the actual Cubit file itself. Hopefully I 
>>> can open it in my older version of Cubit.
>>>
>>> Regards,
>>> J-P
>>>
>>>
>>>
>>> On Thursday, October 20, 2016 at 10:09:47 PM UTC+2, Oded Yaakobi wrote:

 Hello group,

  

 I am interested to test a code that I am developing by solving with it 
 the sp

Re: [deal.II] question on contact algorithm for step-41

[Wolfgang Bangerth]

Hoang,


I have some doubts over the contact algorithm used in step-41. Maybe someone
can help me to clarify:
+ In step-41, instead of solving for delta_U and delta_Lambda, we solve for
U_k+1 by using a finite difference equation. Is it some form of explicit
contact algorithm?


I'm not familiar with the literature on contact problems, so don't know the 
specific terms. You can relate the approach taken in step-41 to what one can 
do when using a the Newton method to minimize an energy functional. There, you 
can solve for the update dx^(k+1) in step k+1 via

  H^k dx^(k+1) = -J^k
(H^k=Hessian, J^k=Jacobian of the energy functional, both linearized around 
the current iterate x^k), but you can also add H^k x^k to both sides to obtain 
the equation

  H^k (x^k+dx^(k+1)) = H^k x^k - J^k
which you can rewrite in terms of the *next* iterate
  H^k x^(k+1) = H^k x^k - J^k

In other words, you can use Newton's method to either solve for the 
*increment* dx^(x+1) or for the *next iterate* x^(x+1). It's still the same 
method.




It is also worth noting that the gap function is not
linearized. In this example, the gap G is constant function by the way,
however considering possible sliding of the node, the gap still need to be
linearized if using the implicit contact algorithm.


I don't quite understand. The gap function is a linear function in this 
example -- it's just the distance between the obstacle and the membrane.




+ The convergence of the contact algorithm is in 18 Newton steps, comparing to
typically 7-8 steps for an implicit contact algorithm. Although it's hard to
compare here because the primal-dual active set is used. But if implicit
contact is used, will the convergence be improved?


I don't know. I think details such as the number of Newton steps required 
critically depend on so many factors (the obstacle function, the energy 
functional, the exact choice of penalty factors, ...) that it's difficult to 
make such comparisons. One would have to try this out by implementing an 
alternative method. If you wanted to try this, we'd be very happy to take your 
program and incorporate it into the tutorial or the code gallery!


Best
 W.

--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Wolfgang Bangerth]

On Thursday, October 20, 2016 at 9:48:55 PM UTC-5, Wolfgang Bangerth wrote:

What is the norm of the first right hand side? You now set the tolerance in

the first iteration to a fixed value of 1e-16, but how does this compare to
the previous value of 1e-12*system_rhs.l2_norm()?


The system_rhs.l2_norm() before the first call to CGSolverat in the zero time
step was 1.82492e-09, so that the 1e-12*system_rhs.l2_norm() became in the
order of 1e-21 by which CGSOlver didn't converge.


OK. I have no suggestion, though, why one or the other would lead to 
convergence or non-convergence of CG.




That still seems wrong to me. As I mentioned in a previous email, you
ought to
make sure that the matrix you build is the matrix you *want* to build.


I was wondering how to make sure the matrix I built is what I meant.


You need to compare properties of your bilinear form with those of your 
matrix. For example, if your bilinear form is symmetric, then your matrix is 
as well -- and that is something you can test. If your bilinear form is 
positive definite, then your matrix needs to be as well and as a consequence 
CG needs to converge. If CG does not seem to converge but your bilinear form 
is s.p.d., you know you have a bug in your matrix assembly. If CG converges on 
one processor but not on two for the same problem, then you know you have a 
bug in the parallel matrix assembly. Etc.


Since you haven't stated your bilinear form or the problem you are solving, I 
can't tell you whether the matrix is supposed to be s.p.d.


Best
 W.


--

Wolfgang Bangerth  email: bange...@colostate.edu
   www: http://www.math.colostate.edu/~bangerth/

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.

To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Re: [deal.II] Re: How to use Trilinos instead of Petsc in step-40

[Hamed Babaei]

Dear Wolfgang,

On Thursday, October 20, 2016 at 9:48:55 PM UTC-5, Wolfgang Bangerth wrote:
>
> What is the norm of the first right hand side? You now set the tolerance 
> in 
>
the first iteration to a fixed value of 1e-16, but how does this compare to 
> the previous value of 1e-12*system_rhs.l2_norm()? 
>

The system_rhs.l2_norm() before the first call to CGSolverat in the zero 
time step was 1.82492e-09, so that the 1e-12*system_rhs.l2_norm() became in 
the order of 1e-21 by which CGSOlver didn't converge.

>
> That still seems wrong to me. As I mentioned in a previous email, you 
> ought to 
> make sure that the matrix you build is the matrix you *want* to build. 
>
 
I was wondering how to make sure the matrix I built is what I meant. 
In addition, Bruno already mentioned sometimes PreconditionerAMG has 
problem with non-positive definite System_matrices. 

Thanks  

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: Error in GridIn::read_abaqus for a mesh with two blocks

[Jean-Paul Pelteret]

Hi Oded,

Two things:
1. The documentation originally stated that "Files generated using Cubit 
11.x, 12.x and 13.x are valid", but I see that its now been updated to 
include 14.x and 15.x. I personally don't know whether this works with 
newer versions of Cubit, but hopefully it has indeed been tested.
2. You've not provided all of the data necessary to recreate your geometry:

 CUBIT> create sphere radius 1
> Successfully created sphere volume 1 
> Journaled Command: create sphere radius 1
> CUBIT> sculpt parallel volume all processors 8 fileroot "sculpt_parallel" 
> gen_input_file "sculpt_parallel.i" overwrite xintervals 8 yintervals 8 
> zintervals 8 box location position -3 -3 -3 location position 3 3 3 smooth 
> 8 num_laplace 2 max_opt_iters 5 opt_threshold 0.6 max_pcol_iters 100 
> pcol_threshold 0.2 sideset 2 void adapt_type 1 adapt_levels 2 combine 
> import show clean
> ERROR: Unrecognized Identifier: 'gen_input_file'


Can you please provide all whatever files are necessary to recreate the 
geometry? Either that or provide the actual Cubit file itself. Hopefully I 
can open it in my older version of Cubit.

Regards,
J-P



On Thursday, October 20, 2016 at 10:09:47 PM UTC+2, Oded Yaakobi wrote:
>
> Hello group,
>
>  
>
> I am interested to test a code that I am developing by solving with it the 
> specific problem of Stokes flow in and around a fluid sphere.   
>
>  
>
> I used Cubit 15.2 to create a mesh of a sphere in a box, as shown in the 
> attached picture, and exported it to an abaqus file. The resulting mesh is 
> composed of two blocks - block 1 for the domain inside the sphere and block 
> 2 for the domain outside the sphere. As a first step, I didn’t define 
> different material id’s to each one of the blocks.
>
>  
>
> When I tried to read the file using GridIn::read_abaqus, I received the 
> following Error message:
>
>  
>
> 
>
> Exception on processing: 
>
>  
>
> 
>
> An error occurred in line <864> of file 
>  in 
> function
>
> void dealii::GridIn<3, 3>::read_abaqus(std::istream &) [dim = 3, 
> spacedim = 3]
>
> The violated condition was: 
>
> false
>
> The name and call sequence of the exception was:
>
> ExcMessage("Internal conversion from ABAQUS file to UCD format was 
> unsuccessful.Are you sure that your 
> ABAQUS mesh file conforms with the 
> requirementslisted in the 
> documentation?")
>
> Additional Information: 
>
> Internal conversion from ABAQUS file to UCD format was 
> unsuccessful.Are you sure that your 
> ABAQUS mesh file conforms with the requirements
> listed in the documentation?
>
> 
>
>  
>
> Aborting!
>
> 
>
>  
>
>  
>
>  
>
> When I tried to save only one of the blocks (either block 1 or block 2) in 
> separate files, there was no problem in reading the corresponding abaqus 
> files.
>
>  
>
> Have I done something wrong in the way I created the multi-block mesh, or 
> is there a bug with GridIn::read_abaqus?
>
>  
>
> Attached are the abaqus files, the Cubit journal file and my code that 
> fails to read the multi-block mesh file.
>
>  
>
> Thank you,
>
> Oded 
>
>  
>

-- 
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see 
https://groups.google.com/d/forum/dealii?hl=en
--- 
You received this message because you are subscribed to the Google Groups 
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email 
to dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[deal.II] Re: How to apply boundary values for a particular point on the boundary instead of the whole boundary surface

[Jean-Paul Pelteret]

Dear Hamed,

The problem here is that VectorTools::interpolate_boundary_values always 
applied boundary values to surfaces in 3d, or lines in 2d. You cannot 
directly apply point constraints using this function. What is happening 
here is that you've set a component mask which says that for ALL of the 
specified components on a given surface, you wish to add some constraint. 
Next you define the constraint using a function. By default, the values 
returned by the function are zero (this is what Vector &values is 
initialised as), and you then can change their value based on position. 
Note that it is the value that changes, not whether they are actively 
constrained or not. So, your function is effectively doing nothing tha the 
ZeroFunction wouldn't already do, since based on some position-dependent 
criterion you're changing the returned values from zero to zero!

However, if you're really wanting to implement point constraints then there 
is some reprise. Look at this post 

 
and other posts 
 
on how to use the vertex_dof_index, and then manually add constraint lines 
your constraint matrix.

Regards,
J-P

On Friday, October 21, 2016 at 1:38:19 AM UTC+2, Hamed Babaei wrote:
>
> Hi friends,
>
> For an elastic problem, I am going to apply zero boundary displacements 
> for three specific points on the center of -x, -y and -z planes of a cubic 
> domain.
> I have already done this but for the boundary surface not a boundary point 
> (the same as incremental_boundary_displacement in step-18). 
> The following is what I wrote to do so which doesn't work properly. In 
> fact, it fixes the whole -x, -y and -z surfaces, not just for the three 
> points on them that I intended. 
>
>   template 
>   class BoundaryCondition :  public Function
>   {
>   public:
>  BoundaryCondition (const int boundary_id);
> virtual void vector_value (const Point &p,
> Vector   &values) const;
> virtual void vector_value_list (const std::vector > 
> &points,
>std::vector >   
> &value_list) const;
>   private:
> const int boundary_id;
>
>   };
>   template 
>   BoundaryCondition::BoundaryCondition (const int boundery_id)
> :
> Function (dim),
> boundary_id(boundary_id)
>
>   {}
>   template 
>   inline
>   void
> BoundaryCondition::vector_value (const Point &p,
> Vector   &values) const
>   {
> Assert (values.size() == dim,
> ExcDimensionMismatch (values.size(), dim));
>
> Point point_x;
> point_x(1) = 5;
> point_x(2) = 5;
>
> Point point_y;
> point_y(0) = 5;
> point_y(2) = 5;
>
> Point point_z;
> point_z(0) = 5;
> point_z(1) = 5;
>
>
> if  (boundary_id ==0 && ((p-point_x).norm_square() 
> <(0.5e-9)*(0.5e-9)))
>values(0) = 0;
> else if (boundary_id ==2 && ((p-point_y).norm_square() < 
> (0.5e-9)*(0.5e-9)))
>values(1)= 0;
> else if (boundary_id ==4 && ((p-point_z).norm_square() < 
> (0.5e-9)*(0.5e-9)))
>values(2)= 0;
>
>   }
>   template 
>   void
> BoundaryCondition::vector_value_list (const 
> std::vector > &points,
>  std::vector > 
>   &value_list) const
>   {
> const unsigned int n_points = points.size();
> Assert (value_list.size() == n_points,
> ExcDimensionMismatch (value_list.size(), n_points));
> for (unsigned int p=0; pBoundaryCondition::vector_value (points[p],
> value_list[p]);
>   }
>
> .and in the constraint I have 
> added VectorTools::interpolate_boundary_values for every boundary plane, 
> for example for the -x plane it looks like :
>
>
> VectorTools::interpolate_boundary_values(dof_handler,
>  boundary_id,
> 
>  BoundaryCondition (boundary_id),
>  constraints,
> 
>  fe.component_mask(x_displacement));
>
>
> I don't know why it doesn't recognize the if condition "  if 
>  (boundary_id ==0 && ((p-point_x).norm_square() <(0.5e-9)*(0.5e-9)))" !!!
>
> I was wondering if you know where I am making mistake, or if there is any 
> step in which this boundary condition has applied.
>
> Thanks,
> Hamed
>

-- 
The deal